Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

DR$\vert$GRADUATE: uncertainty-aware deep learning-based diabetic retinopathy grading in eye fundus images

DR$\vert$GRADUATE: uncertainty-aware deep learning-based diabetic retinopathy grading in eye... Medical Image Analysis (2020) Contents lists available at ScienceDirect Medical Image Analysis journal homepage: www.elsevier.com/locate/media DRjGRADUATE: uncertainty-aware deep learning-based diabetic retinopathy grading in eye fundus images a,b a,b c d,e d d,e Teresa Araujo ´ , Guilherme Aresta , Lu´ ıs Mendonc ¸a , Susana Penas , Carolina Maia , Angela Carneiro , Ana Maria a,b a,b Mendonc ¸a , Aurelio ´ Campilho INESC TEC - Institute for Systems and Computer Engineering, Technology and Science, 4200-465 Porto, Portugal Faculty of Engineering of University of Porto, 4200-465 Porto, Portugal Department of Ophthalmology of Hospital de Braga, Braga, Portugal Department of Ophthalmology of Centro Hospitalar Sao ˜ Joao, ˜ Porto, Portugal Department of Surgery and Physiology of Faculty of Medicine of University of Porto, Porto, Portugal A R T I C L E I N F O A B S T R A C T Diabetic retinopathy (DR) grading is crucial in determining the adequate treatment Article history: and follow up of patients, but the screening process can be tiresome and prone to er- rors. Deep learning approaches have shown promising performance as computer-aided diagnosis (CAD) systems, but their black-box behaviour hinders clinical application. Diabetic retinopathy grading, Deep learning, Uncertainty, Explainability We propose DRjGRADUATE, a novel deep learning-based DR grading CAD system that supports its decision by providing a medically interpretable explanation and an es- timation of how uncertain that prediction is, allowing the ophthalmologist to measure how much that decision should be trusted. We designed DRjGRADUATE taking into account the ordinal nature of the DR grading problem. A novel Gaussian-sampling ap- proach built upon a Multiple Instance Learning framework allow DRjGRADUATE to infer an image grade associated with an explanation map and a prediction uncertainty while being trained only with image-wise labels. DRjGRADUATE was trained on the Kaggle DR detection training set and evaluated across multiple datasets. In DR grading, a quadratic-weighted Cohen’s kappa () between 0.71 and 0.84 was achieved in five dif- ferent datasets. We show that high  values occur for images with low prediction uncer- tainty, thus indicating that this uncertainty is a valid measure of the predictions’ quality. Further, bad quality images are generally associated with higher uncertainties, showing that images not suitable for diagnosis indeed lead to less trustworthy predictions. Ad- ditionally, tests on unfamiliar medical image data types suggest that DRjGRADUATE allows outlier detection. The attention maps generally highlight regions of interest for diagnosis. These results show the great potential of DRjGRADUATE as a second- opinion system in DR severity grading. c 2020 Elsevier B. V. All rights reserved. c 2019. Licensed under the Creative Commons CC-BY-NC-ND 4.0 li- cense http://creativecommons.org/licenses/by-nc-nd/4.0/. guilherme.m.aresta@inesctec.pt (Guilherme Aresta), Corresponding author amendon@fe.up.pt (Ana Maria Mendonc ¸a), campilho@fe.up.pt (Aurelio ´ e-mail: tfaraujo@inesctec.pt (Teresa Araujo), ´ Campilho) arXiv:1910.11777v2 [eess.IV] 29 May 2020 2 Teresa Araujo ´ et al. / Medical Image Analysis (2020) A B C Table 1: International dDiabetic retinopathy (DR) severity scale. NPDR: Non- proliferative DR; PDR: proliferative DR. Non-referable DR: R0 and R1; Refer- able DR: R2, R3 and R4. Microaneurysms Hemorrhages Hard exudates Grade Description R0 - No DR No visible sign of abnormalities Mild R1 - Presence of microaneurysms only NPDR Moderate More than just microaneurysms but Cotton-wool Neovascu- R2 - larization spots NPDR less than severe NPDR Any of the following: Fig. 1: Examples of typical retinal lesions characteristic of diabetic retinopathy. >20 intraretinal hemorrhages Severe R3 - Venous beading NPDR Intraretinal microvascular abnormalities no signs of PDR 1. Introduction Either or both of the following: R4 - PDR Neovascularization Diabetic retinopathy (DR) is a complication of diabetes and Vitreous/pre-retinal hemorrhage is one of the leading causes of blindness worldwide (The Royal College of Ophthalmologists, 2012; Wild, 2004) and the num- ber of diabetic patients is expected to grow from 346 million in 2012 to 552 million people in 2030 (Schully, 2012; The Royal DR is generally classified according to the absence or pres- College of Ophthalmologists, 2012). The majority of visual ence of new abnormal vessels as non-proliferative DR (NPDR) loss cases can be prevented with early detection and adequate or proliferative DR (PDR), respectively. NPDR can be further treatment (Garg and Davis, 2009; Wu et al., 2013). However, graded as mild, moderate and severe (The Royal College of one third of the diabetic patients su er from DR without hav- Ophthalmologists, 2012). These stages can be ordered accord- ing symptoms of visual problems, leading to the progression of ing to their risk of progression (Arev ´ alo et al., 2013). Namely, the disease without treatment (Acharya et al., 2012). Conse- Table 1 represents the International Clinical DR Scale (?), quently, regular check-ups via screening programs are essential adapted to consider only one image for diagnosis (?) by remov- to achieve early diagnosis of DR. ing the quadrant concept, since in the case of a single fundus image analysis quadrant assessment is limited due to the under- DR screening programs have been developed across 1 2 representation of some of the quadrants. Thus, in practice, the the world, namely in the United Kingdom , Ireland exact counting of lesions in each quadrant is not performed by Netherlands (Abramo and Suttorp-Schulten, 2005), 4 5 the specialists, but rather the estimation of what happens in the France (Massin et al., 2008), United States (Cuadros full eye based on their experience and medical knowledge. and Bresnick, 2009), Canada (Nguyen and Wong, 2009) and Portugal (Dutra Medeiros et al., 2015). In these programs, DR Mild NPDR is characterized by the presence of microa- detection and DR severity grading are performed by trained neurysms (MAs), which are the first sign of DR. Hemorrhages specialists by visual analysis of color eye fundus images. (HEMs) start to appear with the evolution of the disease in the These images are captured by fundus photography, allowing deeper layers of the retina. The disease progression results in an for non-invasive visual assessment of several retinal structures, increase of MAs, dot-hemorrhages (red lesions), retinal thick- as depicted in Fig. 1: 1) the fovea, a reddish area without ves- ening, hard exudates (EXs), and the appearance of intraretinal sels located at the center of the macula; 2) the optic disc (OD) microvascular abnormalities (IrMAs), venous beading and cot- which is a round structure close to the macula that corresponds ton-wool spots (CWSs), also known as soft exudates, in the in- to the exit of the axons of ganglion cells that constitute the nermost retina. Finally, PDR results from abnormal vessel for- optic nerve and the door for the retinal vessels that supply the mation, i.e., neovessels (NVs). Other signs of PDR are the vit- retina, and 3) blood vessels, that radiate from the OD center. reous/pre-retinal hemorrhages (PHEMs) and pre-retinal fibrosis During screening, specialists search for abnormalities in these (PFIB) (Hall, 2011). Fig. 1 shows examples of some of these structures and in the retinal tissue and classify the severity of lesions. the disease according to the type and quantity of the findings. Computer-aided diagnosis (CAD) systems can improve the Despite the overall success of these programs, the diagnosis DR screening pipeline both reducing the burden (Scotland et al., process is still prone to errors due to the large number of 2007) and by proving a second objective opinion to the oph- patients, poor image acquisition and variety of lesions. thalmologists, reducing subjectivity in the diagnosis (Gargeya and Leng, 2017; Abramo ` et al., 2016; Gulshan et al., 2016; Quellec et al., 2017). Namely, deep learning has recently al- https://www.gov.uk/topic/population-screening-programmes/ lowed for these systems to achieve near-human performance in diabetic-eye DR detection, i.e. detection of at least one DR-related lesion. https://www.diabeticretinascreen.ie On the other hand, DR grading, i.e. staging of the pathology http://www.eyecheckklant.nl according to its severity, is a more complex problem since it re- reseau-ophdiat.aphp.fr https://www.eyepacs.com quires the identification and integration of di erent abnormal- Teresa Araujo ´ et al. / Medical Image Analysis (2020) 3 ities. In fact, specialists tend to disagree on the grade of com- performed by a 11 convolution and global max-pooling, fol- plex cases (Krause et al., 2018) and thus the external opinion lowing the assumption that the presence of one unhealthy patch of CAD systems may contribute to further improve the success is enough to label the entire image. Because of this, the out- of DR grading. However, due to this higher complexity and the put of this convolution is a map where each element indicates inherent black box behaviour of deep learning systems, it may the malignancy score of a corresponding patch of the input im- be dicult for specialists to understand and accept the system’s age, thus explaining why the image label was predicted. In- decision. With this in mind, we propose a novel DR grading stead, Gonzalez-Gonzalo ´ et al. (2018) has applied an iterative CAD system that not only supports its decision by providing a inpainting approach that, at each step, allows to identify pro- medically interpretable explanation but also estimates how un- gressively less discriminative pixels of the image for DR classi- certain that prediction is. By doing so, the ophthalmologist not fication. This method involves the computation of the saliency only has a reasoning behind the system’s decision, but also a map at each step, i.e., the derivative of the classification output measure of how much that decision should be trusted, making with respect to the input image using guided back-propagation. the DR diagnosis process less prone to errors. The authors use a pre-trained VGG-16 on the Kaggle dataset (512512 pixel images). Similarly to other grading approaches, the problem is treated as being ordinal and thus the model is op- 1.1. State-of-the-art timized using a mean-square error loss. To avoid overfitting due Deep learning has become the preferred approach, over field to the high class imbalance, the classes are balanced taking into knowledge-based feature engineering approaches, for DR grad- account the referal/non-referal sample distribution. ing since it allows to better exploit the large number of data Besides explaining the decision, uncertainty estimation of available and to better deal with the labeling noise resulting a model’s prediction is of interest as it helps to avoid conse- from the complexity of the task. For instance, Krause et al. quences of the blind use of the model response (Kendall and (2018) fine-tuned an Inception-v4 with 1,600,000 images of Gal, 2017). This is particularly true on medical settings, where size 779779 pixels to avoid the elimination of small lesions misdiagnosis can have severe consequences for the patients. such as MAs. Also, to account for data imbalance, the authors With this in mind, Leibig et al. (2017) analysed the uncer- apply a cascade of thresholds on the output probabilities of the tainty information from deep neural nets for DR detection and model for each DR grade to obtain the final image classification. referable DR detection. The authors tested the dropout-based These thresholds are sequentially applied in a severity descend- Bayesian uncertainty estimation, comparing it with alternative ing order to ensure that a minimum grade-wise sensitivity is techniques, such as directly analysing the network’s softmax achieved. Instead, de la Torre et al. (2018) considered DR grad- output. In the dropout-based approach, dropout is switched ing as an ordinal classification problem and thus aimed at reduc- on during inference. Performing multiple predictions on the ing the distance between predicted and true labels of an image. same input image can be interpreted as having an ensemble of For that, they proposed a quadratic-weighted-Kappa-based loss highly similar models, and thus disagreement between the in- function, which allowed to improve the performance of a DR ferences is indicative of the global model uncertainty (Gal and grading network comparing to training with the cross-entropy Ghahramani, 2016). The authors state that the Bayesian treat- loss. The authors tested di erent input images sizes (squares ment worked better than the softmax output as a proxi for un- with side of 128, 256, 384 and 512 pixels), and the highest res- certainty and also show that uncertainty-aware decision referral olution yielded the best performance. The network, which has can improve the diagnosis. 11.3 million parameters for the 512512 case, was trained in a public dataset (Kaggle DR detection). To overcome the high unbalancing of the dataset, the data is artificially balanced using 1.2. Contributions data augmentation, oversampling the least represented classes. Despite the high classification performance of these systems, Many research work has already addressed DR grading, ex- their black box behavior hinders the application on screening plainability and uncertainty estimation separately. Instead, in settings. Indeed, providing explanations for DR grading is chal- this study we propose DRjGRADUATE, a system that is capa- lenging due to the complexity of the task, and thus research ble of simultaneously providing an explanation and a measure is mainly focused on the binary detection and referral tasks. of uncertainty together with a DR severity grade label. Our For instance, Quellec et al. (2017) analyzed several explana- main contributions are the following: tion methods for the referable DR decisions of their CNNs’. Namely, the authors compare deconvolution (Zeiler and Fergus, estimation of a DR grade with an associated uncertainty, 2014), sensitivity analysis (Simonyan et al., 2014) and layer- thus allowing to better establish which cases require fur- wise relevance propagation Bach et al. (2015). Based on (Si- ther inspection by specialists, which is of high interest in monyan et al., 2014), the authors also proposed the creation of computer-aided diagnosis systems; a pixel-wise explanation map that is jointly optimized with the CNN prediction during training. Costa et al. (2019) proposed  explanation of the algorithm decision, easily and directly the attention map generation via a pre-trained network coupled interpretable by the retinal specialists, in the form of an at- with a Multiple Instance Learning (MIL)-based approach for tention map, which will leverage its use by medical profes- DR detection. Namely, the network is cropped at early layers sionals since it partially mitigates the black-box behaviour to reduce the model’s receptive field, and the classification is associated with deep nets; 4 Teresa Araujo ´ et al. / Medical Image Analysis (2020) extensively validated deep learning-based method for DR 2.1. Predicting a grade with uncertainty severity grading of eye fundus images in the internation- For each input image I, DRjGRADUATE predicts a DR ally recognized DR levels. grade y ˆ and a grade uncertainty u. Let M be the output of the last layer of DRjGRADUATE’s backbone after assessing an input image I. M has size N  N  F, where N and F are, 2. DRjGRADUATE for DR grading respectively, the side and number of features. The lesion map L is computed by performing a 1  1  1 convolution with a The backbone of DRjGRADUATE, depicted in Fig. 2, is linear activation over M, so that L 2 R has size N  N  1. composed of several convolutional-batch normalization blocks L is a map where the value of each element L indicates the i j interleaved with max-pooling layers that has already shown po- y presence of a lesion of grade bL e on a patch of size depen- i j tential for DR grading tasks (de la Torre et al., 2017). The top dent on the model’s receptive field in the input. Following the of the network is composed of a novel component that allows MIL approach, the grade of an image is initially estimated as for the inference of the image’s grade y ˆ together with the cor- y ˆ = max(L) (Costa et al., 2019). responding explainabilityE and the explicit uncertainty u of the Having into account assumptions 1 and 3, y ˆ 7! y ˆ must r g decision, without using labels other than the DR grade, such as ensure that p > p > p > : : : , and p > p > by ˆ c by ˆ c1 by ˆ c2 dy ˆ e dy ˆ e+1 r r r r r lesion-wise annotations. p > : : : , i.e., that the classes near to y ˆ have higher prob- dy ˆ e+2 r Let DRjGRADUATE be a classification model that predicts ability than the others. For each image, DRjGRADUATE must y ˆ = argmax (p), where p is the class-wise probability and g c thus predict a generalized Bernoulli distribution biased around p = 1, with C denoting the number of considered classes. c the classes closer to y ˆ . This bias is introduced in the model by The following assumptions regarding DR grading are made: means of a Gaussian distribution N centered on y ˆ and thus the [assumption 1] grading is a discrete ordinal classification image-wise grade probability p is: problem with an image label y 2 G = f0; 1; 2; 3; 4g, where G is the set of possible DR grades; in this scenario, we want ( ! ) 1 [i y ˆ ] to optimize DRjGRADUATE by minimizing the distance be- r p = exp 0:5 ; 8c 2 G tween the target and the prediction, min(jy y ˆ j), and [assump g 2 tion 2] there is always a di erent lesion type per grade on the p = P (1) universe of possible DR-related lesions, DR ; thus, if one lesions c2G c associates a numerical value to each lesion from DR , ac- lesions cordingly to the DR stage in which they appear, and thus y ˆ where  is the variance of N and p 2 [0 ; 1]. Fig. 3a illus- can be predicted by following the standard MIL assumption, trates the sampling of the probability values p from the Gaus- i.e. max (DR ) 7! y ˆ . lesions g sian probability distribution. Knowing that the uncertainty of DR grading can be seen as a regression problem and con- the prediction can be measured by the value of its entropy (Fris- sequently the predicted image grade y ˆ can be inferred from a ton, 2010), i.e. u is directly related to S = p log(p ), it fol- c c 2 2 continuous output y ˆ 2 R . This also allows for the implicit r + lows that since S depends on  , then u depends on  . Thus Multiple Instance Learning (MIL) of the problem, as we can to by estimating  it is possible to infer the uncertainty associated retrieve y ˆ from a predicted lesion map L, y ˆ = max(L). Further- r r with each predicted grade, y ˆ = argmax(p). more, the uncertainty u could be estimated by u / jby ˆ e y ˆ j . r r In DRjGRADUATE,  is inferred by using a fully- However, this uncertainty does not properly hold in scenarios connected layer on top of M. Since both y ˆ and  de- where the regression cut-o thresholds are adjusted to better pend on M, both grading and the associated uncertainty af- fit the data distribution and malignancy priority, as in many fect the backbone of the model, easing its convergence. In- Kaggle competition solutions . Also, in terms of clinical ap- deed, DRjGRADUATE is trained end-to-end minimizing the plicability, it is of higher interest to provide a class-wise prob- loss function: ability to clinicians since it is easier to interpret. Because of this, [assumption 3] DRjGRADUATE assumes that the pre- L = [b log(p )] + (1 ) (2) c c diction and its uncertainty can be modeled by a Gaussian curve c2G centered near the most probable class. Namely, for each im- where b is the hot-encoded (i.e., a binary representation of) i age, DRjGRADUATE predicts a variance  , which allows for and 2 [0; 1] is a weighting factor. The first term, the categori- y ˆ 2 R 7! y ˆ = argmax(p). With this, DRjGRADUATE is r + g cal cross-entropy, leads DRjGRADUATE to predict y ˆ and thus 2 r capable of providing an uncertainty  and a class-wise prob- y ˆ as close to the truth as possible. The second term regularizes ability that still attempts at minimizing the distance between to avoid excessive Gaussian spreads and thus the stagnation the prediction and the ground truth. The next sections describe of the learning on the local minimum of near equiprobable pre- DRjGRADUATE with more detail. dictions. The global loss L and the structure of DRjGRADUATE im- plicitly introduce an inter-dependency between y ˆ and  that be denotes rounding to the nearest integer, bc rounding down and de allows to deal with images of di erent complexity. For less rounding up. complex images, both  and the preliminary class prediction https://storage.googleapis.com/kaggle-forum-message- attachments/88655/2795/competitionreport.pdf error, jby ˆ e y ˆ j, should be low. However, for dicult images, r r Teresa Araujo ´ et al. / Medical Image Analysis (2020) 5 MP Lesions Conv ... 2×(1×1) N N 1 640 64 FC Uncertainty Maxpool 2×2 Block i Convolution 3×3 + Batch norm + ReLU Convolution 3×3 + Batch norm + ReLU Fig. 2: DRjGRADUATE’s architecture. The network is composed of several convolutional-batch normalization blocks interleaved with max-pooling layers. For each input image I, DRjGRADUATE predicts a diabetic retinopathy (DR) grade y ˆ and a grade uncertainty u. M is the output of the last layer of DRjGRADUATE’s backbone after assessing I. The lesion map L is computed by performing a 1  1  1 convolution with a linear activation over M. L is a map where the value of each element L indicates the presence of a lesion of grade bL e on a patch of size dependent on the model’s receptive field in the input. Following the Multiple i j i j Instance Learning approach, the grade of an image is initially estimated as y ˆ = max(L). The output probabilities (p ; c 2 f0:::4g) are then computed from y ˆ based r c r on a Gaussian distribution centered on y ˆ with variance equal to the uncertainty u, with u learned during training. high values of  can compensate for classification errors since two-dimensional multivariate normal distribution with standard 0 0 the higher spread of the Gaussian tends to increase the predicted deviation  , with  related to the receptive field of the net- probability of the correct class, lowering the cross-entropy loss. work, and s is the stride of DRjGRADUATE. The convolution Fig. 3b shows the evolution of L (Eq. 2) as function of y ˆ and with the Gaussian kernel allows to give higher priority to central , for di erent targets y, y 2 G. These plots show that when pixels of the receptive field, which have a higher contribution to y ˆ is close to y the L loss reaches its minimum. Further, if one the value of each L . r i j freezes y ˆ at a value close to the ground truth, the  that leads 2.3. Dealing with class imbalance to a lowerL loss is smaller than if y ˆ is fixed at a value far from DR grading is an imbalanced problem skewed towards the the ground truth. healthy cases, similarly to many other medical imaging prob- lems. Data imbalance introduces convergence issues during 2.2. Explainability of the DR grade model training, demanding for class balancing strategies. Over- The goal of this study does not include a pixel-level pre- sampling has shown to be an overall good choice for dealing diction of lesions. Even human-performed lesion segmenta- with class imbalance (Buda et al., 2018). However, exces- tion shows high variability, due to the extremely small size of sively reducing class imbalance may lead to an unnecessarily some lesions and also to di erences in the annotation style. In- high trade-o between the model’s sensibility and specificity, deed, some ophthalmologists prefer to only annotate clusters i.e. increase the misclassifications of healthy images, since the of lesions, whereas others are more detailed and prefer to give most represented classes will now tend to be relatively over- lesion-wise markings. Since DR image grading is performed classified. With this in mind, DRjGRADUATE is trained using based on the presence of the typical DR lesions, pinpointing an iterative batch balancing scheme that has already shown suc- relevant lesions should be a valuable tool for doctors to assess cess for DR grading : the algorithm’s prediction without the need for a refined seg- mentation. t1 t1 w = r w + (1 r )w (4) ct co c f DRjGRADUATE provides a grade-wise explanation map, as where w is the expected ratio of images of grade c on the batch ct illustrated in Fig. 4. For each grade c the lesion presence en- at epoch t, w is the representativeness of class c on the dataset, c 0 coded in L is translated to an explainability map E of the same w is the ratio of class c on the batch after f training epochs c f size of I via: and r is the decay rate that controls the update of w . With this, ct in the beginning of the training the network is confronted with > a completely balanced set and does the initial learning based on 1 if c 0:5  L < c + 0:5 i j E (s i; s j) = this, and then progressively starts to adjust (or fine-tune) to a 0 otherwise more real data distribution, but still not as skewed as the origi- E = E ~N 0 (3) c 2 ; nal training set distribution. The oversampling of the least rep- resented classes is performed via online horizontal and vertical where E is a matrix of the same size of E , c 2 G n f0g is one of the possible pathological grades, i; j are each row and column position, respectively, in the L map of side N https://www.kaggle.com/blobs/download/ (i; j 2 N j i; j < N), and N is a convolutional kernel from a forum-message-attachment-files/2797/report.pdf 2 ; Block 1 Block 5 6 Teresa Araujo ´ et al. / Medical Image Analysis (2020) y = 0 y = 1 y = 2 y = 3 y = 4 0.0 0.16 0.31 p 0.47 0.63 0 1 ŷ 2 3 4 Diabetic retinopathy grade 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 (a) Probability Gaussian distribution with ŷ ŷ ŷ ŷ ŷ r r r r r mean y ˆ and standard deviation  (Eq. 1). (b) L loss values (Eq. 2, with = 0:7) as function of y ˆ and , for y 2 G. Loss colorbar: 0 2 . Fig. 3: Uncertainty estimation (Eq. 1) and loss representation (Eq. 2). Messidor-2 Kaggle (test) IDRID SCREEN-DR DMR 0.5<L ≤1.5 Gaussian processing 1.5<L ≤2.5 Gaussian processing N Fig. 5: Class distribution of the DR grading datasets used for evaluation. R0/no DR, R1, R2,  R3, R4/PDR,  SDR,  PPDR. 2.5<L ≤3.5 Gaussian N processing Lesions have been used for evaluation. Fig. 5 shows the class distribu- L>3.5 Gaussian processing tion of these datasets. Kaggle DR detection: This dataset was published by Kag- gle , comprising a large number of high resolution images (ap- Fig. 4: Attention maps construction. Gaussian processing is described in Eq. 3. proximately 35 000 in the training set, 55 000 in the test set). Images were acquired under a large variety of conditions, using di erent types of cameras. Images are labeled by a single clin- flips, rotations, brightness adjustments and contrast normaliza- ician with the respective DR grade, out of 4 severity levels: 1 tion. - mild, 2 - moderate, 3 - severe, and 4 - proliferate DR. Level 0 stands for the absence of DR. The organizers state that this 3. Experiments dataset is similar to any real-world dataset, in which noise is present in both images and labels. Consequently, images may Intra- and inter-datasets experiments were performed for contain artifacts, be out of focus and/or with inadequate expo- evaluating DRjGRADUATE in terms of: 1) DR grading per- sition, aiming at promoting the development of robust and re- formance, 2) uncertainty estimation, and 3) explainability. The liable algorithms that behave well in the presence of noise and DR grading performance was quantitatively evaluated on sev- across a variety of images. eral DR-labeled eye fundus image datasets, and class-wise per- Messidor-2: The Messidor-2 public dataset (Abramo ` et al., formance was analysed. The relation between the uncertainty 2016) consists of color eye fundus images, one fovea-centered and the algorithm’s performance was assessed. Further, the un- image per eye, of 874 subjects with diabetes, for a total of 1748 certainty estimation on good and bad quality images was com- images. Images were acquired using a color video 3CCD cam- pared, using image quality-annotated datasets, in order to assess era on a Topcon TRC NW6 non-mydriatic fundus camera with a whether bad quality images would be associated with higher un- 45 FOV, at 1440960, 22401488, or 23041536 pixels (px). certainties. Finally, uncertainty estimation on unfamiliar data The ground truth (Krause et al., 2018) corresponds to an ad- types was also analysed, aiming at inspecting whether the es- judicated consensus of three specialists. A total of 1744 images timated uncertainties could allow to detect unknown image were used, as the remaining were adjudicated as ungradable. types (i.e., outliers). Di erent medical images were consid- IDRID: The Indian Diabetic Retinopathy Image Dataset ered: colon (colonoscopy), skin, cataracts surgery and breast (IDRiD) dataset (Porwal et al., 2018) is composed of 413 im- biopsy microscopy images. For the explainability, the attention ages of 42882848 px, 50 FOV, acquired using a Kowa VX-10 maps produced by DRjGRADUATE were qualitatively evalu- alpha digital fundus camera. Images were classified accord- ated through visual inspection, as well as quantitatively, by ingly to the International Clinical Diabetic Retinopathy Scale. comparison with lesion annotations performed by an ophthal- mologist. 3.1. Datasets https://www.kaggle.com/c/diabetic-retinopathy-detection https://medicine.uiowa.edu/eye/abramo http://www.adcis.net/en/third- 3.1.1. DR grading datasets party/messidor2/ The Kaggle DR detection training dataset was used for train- https://www.kaggle.com/google-brain/ ing DRjGRADUATE, and several public and private datasets messidor2-dr-grades Probability Teresa Araujo ´ et al. / Medical Image Analysis (2020) 7 Table 2: Modified Davis grading used in Takahashi et al. (2017) for the DMR dataset annotation. Grade Description No DR No signs of diabetic retinopathy Microaneurysm, retinal hemorrhage, hard exudate, SDR retinal edema, and more than 3 small soft exudates Soft exudate, varicose veins, intraretinal microvascular PPDR abnormality, and non-perfusion area over one disc area Neovascularization, pre-retinal hemorrhage, vitreous PDR hemorrhage, fibrovascular proliferative membrane, and tractional retinal detachment Fig. 6: Examples of ophthalmologist’s annotations on the SCREEN-DR dataset.  MA, HEM,  EX, CWS, IrMA,  NV,  PHEM,  PFIB. DMR: The dataset announced in Takahashi et al. (2017) (DMR) is composed of 9939 posterior pole color fundus images DR1: This dataset (Pires et al., 2012) is constituted by 5776 (27202720 px) from 2740 diabetic patients. Images were cap- images with resolution 640480 px, captured using a Topcon tured using a NIDEK AFC-30 fundus camera, and have a 45 TRC-50X mydriatic camera (45 FOV). From these, 1300 are FOV. One to four images may exist from a given patient. Im- labeled as being of good quality (do not contain blur and are ages were graded by modified Davis grading, described in Table centered on the macula) and 1392 have poor quality (blur). 2. The international scale grades (Table 1) were converted to the DRIMDB: The Diabetic Retinopathy Image Database modified Davis scale (Table 2) using the following assumption: (DRIMDB) (Sevik et al., 2014) is composed of 216 images di- NDR = R0, SDR = R1 [ R2, PPDR = R3 and PDR = R4. vided into three classes: good (125), poor (69), and outlier (22). SCREEN-DR: This private dataset consists of retinal images The images were captured by a Canon CF-60UVi fundus cam- from a portuguese DR screening program, managed by the Por- era (FOV of 60 ), and have a resolution of 570760 px. Images tuguese North Health Administration (ARSN). A subset of 966 were annotated by an expert, and good-quality images are suit- images was selected (SCREEN-DR dataset). These were ac- able for use in medical diagnosis by an ophthalmologist. quired using Canon CR-2 AF nonmydriatic retinal cameras and HRF: The High-Resolution Fundus (HRF) database (Kohler ¨ have di erent resolutions (with width and height ranging ap- et al., 2013) contains 18 pairs of images from the same eye proximately from 1500 px to 2600 px). Images from both left from di erent subjects, captured using a Canon CR-1 fundus and right eyes were included, which may be centered at the camera (FOV of 45 ). In each pair there is a poor and a good macula, at the optic disc or elsewhere. Images were graded in quality image. These poor quality images have decreased sharp- 5 levels by a retinal specialist. For 348 of these images, pixel- ness, which can be local or global, due, for instance, to a defo- wise annotations of the lesions MAs, HEMs, CWSs, IRMAs, cused camera. Other quality features, such as image contrast or EXs, NVs, PHEMs and PFIB are available. MAs appear as red- illumination conditions, are not contemplated. dish small and circular dots (Zhang et al., 2010), and may oc- cur alone or in groups (Mahendran and Dhanasekaran, 2015). 3.2. Evaluation HEMs are bright red spots and, contrarily to MAs, their shape 3.2.1. DR grading and appearance has substantial variability (Sil Kar and Maity, 2017). Exudates present a large variety of shapes, sizes and The DR grading performance of DRjGRADUATE is eval- contrast (Harangi and Hajdu, 2013). EXs have a more yellow- uated using the Cohen’s quadratic weighted Kappa (), used ish and bright appearance and sharp edges, whereas CWSs are for measuring inter-rating agreement between raters in ordinal more pale yellow/white without well defined borders (Prentasic multi-class problems. This metric penalizes discrepancies be- and Loncaric, 2016). NVs are thin and weak, exhibiting irregu- tween ratings, in a manner that depends quadratically on the lar appearance such as tortuous patterns and loops (Gupta et al., distance between the prediction and the ground-truth, as fol- lows: 2017). Di erent tools were available for annotating the images, P P C C w O i; j i; j such as dots, circles, polygons and free-hand tool, as exempli- i j = 1 (5) P P C C fied in Fig. 6. The annotated lesions were grouped in the di er- w E i; j i; j i j ent grades to allow a direct correspondence with the predicted where C is the number of classes, w is the weight matrix, attention map for each grade. The R1 ground truth map in- O is the observed matrix and E the expected matrix. w is i; j cludes only the annotated MAs; R2 contains only HEMs; R3 the weight penalization for i; j element, which is given by: contains HEMs, CWSs and IRMAs; R4 includes NVs, PHEMs (i j) th . O is the number of samples from the j class (ground and PFIB. EXs were also annotated but not included in the 2 i; j (C1) th grade maps since they are not directly evaluated in the grading truth) which are classified by the model as being from the i process. class. E is the outer product between the prediction and i; j ground truth classification histogram vectors, normalized such 3.1.2. Quality assessment datasets Three di erent eye fundus images datasets with image qual- ity information were used on uncertainty-related experiments. https://www5.cs.fau.de/research/data/fundus-images/ 8 Teresa Araujo ´ et al. / Medical Image Analysis (2020) P P that E = O .  values range from 1 (complete dis- intensity to a progressively lower value, and evaluating the in- i; j i; j i; j i; j agreement) to 1 (perfect agreement). fluence on the prediction’s uncertainty. However,  alone does not completely reflect the perfor- mance of a DR grading algorithm. When dealing with a very 3.2.3. Explainability unbalanced problem, as DR grading (Fig. 5),  is dominated by With the purpose of evaluating the explanation map, a pre- the most represented classes. Using confusion matrices allows dicted object e corresponds to each connected component from to perform grade-wise evaluation and thus understand better the the grade-wise explanation E after applying a fixed probabil- algorithm’s behaviour. We computed normalized versions of ity threshold. The explanation was quantitatively evaluated as the confusion matrices by dividing each value of the original follows: 1) percentage of all e that overlap with the ground matrix by the sum of the values of each row. truth map of the corresponding grade (O ); 2) percentage of ob j g all e that overlap with the ground truth map of the correspond- MIL assumption validation. Preliminary experiments were ing grade or lower (O ); 3) percentage of correctly predicted ob j performed to test whether the MIL assumption hampered the images for which the maximum activation object overlaps with DR grading performance of the network. This was done by the ground truth map of the corresponding grade (O ); 4) per- max training a network in which the MP layer of DRjGRADUATE centage of correctly predicted images for which at least one e (Fig. 2) is replaced by a fully connected (FC) layer, keeping overlaps with the ground truth map of the corresponding grade the remaining network as it was. No significant performance (O ); 5) percentage of ground truth objects that overlap with class difference was found between the non-MIL and MIL models, the predicted map of the corresponding grade or higher (O ); gt with the first tending to present better R0 detection but a higher 6) percentage of all e that overlap with the ground truth map confusion among the disease cases. Because of this, no further of any grade. (O ). Note that the evaluations 2) and 5) are any experiments for this model were performed. justified by the max-pooling assumption made in our pipeline: it is acceptable to assume that if a given region is detected as being from grade 2, for instance, it may also contain signs from 3.2.2. Uncertainty estimation a lower grade which are not highlighted since only the higher Experiments were performed in order to assess how the pre- grade prevails after the max-pooling. diction performance relates with the uncertainty estimation of those predictions. In order to test whether low uncertainty 3.2.4. Statistical tests images are indeed associated with a higher performance, we defined a set of thresholds on the uncertainty values and, for Statistical tests were performed to assess significant di er- each threshold, select the subset of images with estimated un- ences between the uncertainties inferred for di erent datasets. certainty below that threshold. We then compute the  of our Since these data samples did not follow a normal distribu- model predictions on these subsets. Since very low thresholds tion, which is an assumption on most parametric statistical hy- may lead to the inclusion of a number of images which do not pothesis tests, the non-parametric Kruskal-Wallis H Test was truly represent the dataset, and the  computation for these few used (Kruskal et al., 2012). Further analysing the p-value may samples would not be reliable, the uncertainty thresholds start at not be enough when one is dealing with suciently large sam- 0.15 to ensure an acceptable image subset size. This threshold ples, since a statistical test will indicate a significant di erence was determined considering the dataset with the smaller num- in most of the cases (Sullivan and Feinn, 2012), i.e. very large ber of images (SCREEN-DR), for which an uncertainty thresh- sample sizes may lead to the detection of di erences that are old of 0.15 lead to the inclusion of 5% of the dataset, corre- quite small and possibly non relevant. It is thus important to re- sponding to 48 images. Furthermore, cumulative histograms of port the p-value results but also the e ect size, which measures the class-wise number of images as function of the uncertainty the magnitude of the di erence between groups. The e ect size were produced. This allows to exclude the possibility of varia- is defined as ES =   , where  and  are the averages 1 2 1 2 tions in  being influenced by the different proportion of images of the two groups. The Cohen’s d was used as an e ect size per grade instead of being related with the uncertainty. measure, and can be defined as follows: Finally, matrices of the estimated uncertainty were computed by averaging the DRjGRADUATE’s produced uncertainties for d = ES=S (6) all images from each position of the confusion matrix (predicted where S , the pooled standard deviation, is given by: grade vs ground truth grade). 2 2 (n 1) s + (n 1) s Sensitivity analysis. A sensitivity analysis was performed to 1 2 1 2 S = (7) assess the image quality’s influence in the uncertainty esti- n + n 2 1 2 mation. Namely, the influence of image blur was studied by where n and n are the sample sizes for the two groups and s smoothing the input image with Gaussian filters of increas- 1 2 1 and s are the standard deviations for the two groups. ing standard deviation values and inspecting the effect on DRjGRADUATE’s uncertainty estimation. Likewise, contrast The interpretation of d, considering the thresholds proposed was studied since lower contrasts are also characteristic of bad by Cohen and revised by Sawilowsky (2009), can be established quality images. For that, the image’s contrast was decreased via as follows: 1) d = 0:1: very small 2) d = 0:2: small 3) d = 0:5: linear histogram stretching, by setting the image’s maximum medium 4) d = 0:8: large 5) d = 1:3: very large 6) d = 2: huge. Teresa Araujo ´ et al. / Medical Image Analysis (2020) 9 Having into account that our datasets are commonly on the 4.2. Uncertainty estimation order of more than 1000 samples, the e ect size is reported to- The evolution of  with the uncertainty level threshold for gether with the significance level. Since the significant di er- four di erent datasets is shown in Fig. 8a. For the Kaggle DR ence will most likely be achieved in these datasets, d provides detection test set, a  of approximately 0.82 is achieved for the a measure of the magnitude of these di erences. 50% lowest uncertainty images. Likewise, considering the 75% lowest uncertainty predictions,  is of approximately 0.8. The average uncertainty matrix is shown in Fig. 8c, Fig. 8b depicts 3.3. Parameter setting the grade-wise normalized cumulative histogram as function of the inferred uncertainty. Finally, Fig. 9 shows examples of Kag- The input size to the network (Fig. 2) is 640640 3, and the gle DR detection test set images along with DRjGRADUATE’s side of the output feature map, N, is 20. The receptive field, RF uncertainties of the predictions. of the network is of 156, and the stride, s, is of 32. Regarding 0 The uncertainty estimation results for three quality-labeled the attention maps generation, the  of N 0 from Eq.3 was 2; datasets (section 3.1.2) are shown in Fig. 10a. For all datasets, set to RF/6, since it allows to fit the majority of the Gaussian the inferred uncertainty is statistically higher for the bad quality distribution in the kernel window. The weight in the loss subsets (p-value < 0:05), and the respective Cohen’s d is huge (Eq.2) was set to 0.7, in order to preserve the greater importance for the DRIMDB and medium/large for the other datasets. of the log-loss part but still impose some level of regularization. Fig. 10b shows the uncertainty behavior of DRjGRADUATE Regarding the class balancing, we chose r = 0.99 and w = for datasets of unfamiliar data types along with the eye fun- (0:5; 2; 2; 3; 3) ( f = 300) in Eq. 4, aiming at not overfitting the dus image datasets used for evaluation. Between the unfamiliar original dataset distribution and to increase the generalization data types the statistical hypothesis tests indicate that, for each capability of the network. dataset, the results are statistically di erent (p-value < 0:05) from the neighbouring dataset, and the Cohen’s d values sug- 3.4. Training details gest a large to very large magnitude of the di erences. From the BACH to the SCREEN-DR dataset a statistical di erence is Images were cropped around the FOV and resized to the found, with Cohen’s d indicating a medium magnitude of the input size of the network, 640 640 pixels. The model was di erences. Regarding the eye fundus images, some datasets trained for 240 epochs with the Kaggle DR detection training show a statistically di erent average uncertainty. However, in set, using a batch size of 30, the Adam optimizer, and a learning these cases, Cohen’s d values indicate a small to very small rate of 310 . Experiments were performed on an Intel Core magnitude of the di erences. i7-5960X, 32Gb RAM, 2GTX1080 desktop with Python 3.5, Keras 2.2 and TensorFlow 1.8. 4.2.1. Sensitivity analysis Fig. 11 shows examples of images blurred with the Gaus- sian filters of increasing standard deviation (std). Fig. 12 de- 4. Results picts the results obtained for this test in 5 different datasets (Kaggle DR detection test set, Messidor-2, IDRID, SCREEN- DR and DMR), where the curve represents the average of the 4.1. DR grading uncertainty values for each dataset, and the shaded area corre- sponds to a confidence interval of 95%. The sensitivity analysis DRjGRADUATE achieves similar  to other state-of-the-art via contrast adjustment did not suggest any relation between methods, as depicted in Table 3. Note that the majority of image contrast and uncertainty estimation. the state-of-the-art methods did not test on di erent datasets, and thus an assessment of their generalization capabilities is 4.3. Explainability not possible. DRjGRADUATE was evaluated in di erent test datasets and the consistency of the results suggests that the sys- Several examples of DRjGRADUATE’s attention maps vs tem has properly learned discriminant features for DR grading. the specialist’s annotation on SCREEN-DR dataset images are The grade-wise normalized confusion matrices (in percentage) shown in Fig. 13. In each predicted map, the region with max- are shown in Fig. 7. In general, all the matrices have a diagonal imum activation is surrounded by a square. Figs. 13a, 13c, tendency, i.e., the predictions are close to the ground truth. R4 is 13e and 13f illustrate examples where the maximum value of the most misclassified grade. For instance, for both the IDRID the explanation map properly explains the predicted grade. The and SCREEN-DR datasets (Fig. 7c and 7d), the majority of the remaining images show cases of incorrect predictions and/or R4 images is confounded with R3. From the misclassifications incorrect explanations. reported in Fig. 7, the percentage of over-diagnosis is between Table 4 shows the explanation-wise quantitative performance 24% (Kaggle DR detection test set) and 55% (SCREEN-DR) in terms of the number and type of highlighted lesions, O , ob j g for all datasets, which is higher than the 21% expected for hu- O , O , O , O and O (refer to Section 3.2.3). Fig. 14 ob j max class gt any man observers (Krause et al., 2018). Consequently, the per- shows the values of O and O obtained through the appli- max class centage of under-diagnosis is between 45% (SCREEN-DR) and cation of di erent thresholds to the attention maps. The higher 76% (Kaggle DR detection test set), being always lower than the threshold the smaller the objects on the map, resulting in the 79% reported for human observers. higher specificity and lower sensitivity. 10 Teresa Araujo ´ et al. / Medical Image Analysis (2020) Table 3: Quadratic weighted kappa for DRjGRADUATE’s DR grading in four di erent datasets. / 7025 images from the Kaggle DR detection training dataset, * 7000 images from the Kaggle DR detection training dataset, ? 4 classes,  496 out of the 9443 images. Kaggle test Kaggle Messidor-2 IDRID SCREEN-DR DMR? DRjGRADUATE (proposed) 0.74 0.75/ 0.71 0.84 0.74 0.78 de la Torre et al. (2018) 0.74 - - - - - Gonzalez-Gonzalo ´ et al. (2018) - 0.72* - - - - Takahashi et al. (2017) - - - - - 0.79 Predicted Predicted Predicted Predicted Predicted R0 R1 R2 R3 R4 R0 R1 R2 R3 R4 R0 R1 R2 R3 R4 R0 R1 R2 R3 R4 R0 R1 R2 R3 54 40 81 16 3 0 0 4 1 1 79 16 3 2 0 70 21 7 2 0 78 21 1 0 35 48 17 0 0 16 66 18 0 0 25 40 30 5 0 9 41 41 8 1 22 71 6 1 11 25 49 14 1 2 18 62 18 0 0 6 43 49 2 1 6 45 42 6 5 52 38 5 34 55 4 21 73 6 0 0 10 70 20 0 0 15 72 13 0 14 18 68 2 5 0 0 0 17 40 43 0 0 12 59 29 1 5 30 47 17 2 8 22 33 35 0 (a) Kaggle DR detection test (b) Messidor-2 dataset (c) IDRID dataset (d) SCREEN-DR dataset (e) DMR dataset Fig. 7: Confusion matrixes of DRjGRADUATE’s predictions vs the ground truth DR grades for di erent test datasets. Matrixes were normalized by dividing each value of the original confusion matrix by the sum of the values from the corresponding row (total number of images from that grade in the dataset), and the shown values correspond to percentages. Table 4: Assessment of the produced attention maps (threshold=0.3) in com- Because of this, when inferring images with di erent annota- parison with specialist annotated lesions. These metrics are described in sec- tion criteria or with the presence of artifacts, DRjGRADUATE tion 3.2.3. misclassifies images as can be observed in Fig. 15c and 15d. O O O O O O In other scenarios, misclassifications occur due to failed le- ob j g ob j max class gt any sion detection or similarities between di erent lesion types. 0.506 0.677 0.712 0.784 0.526 0.729 Namely, failures such as predicting a R1 image as R0 can be due to missing a single MA. Indeed, the small size and sub- The threshold 0.3 seems to be a good compromise since it al- tlety of these lesions make their detection particularly chal- lows to retain a large percentage of the lesions while producing lenging (Fig. 16a, where no apparent lesion is seen) and easily confoundable with acquisition or anatomical noise (Fig. 16b). regions not too coarse. Confusions between R1 and R2 are likely related with the di- culty of distinguishing between the red lesions MAs and HEMs 5. Discussion (Fig. 16c and 16d), a task that also challenges medical experts since they can be very similar in size and color. Furthermore, 5.1. DR grading di erent grading standards also account for these classifica- tion errors. Namely, for the DMR dataset, the conversion of DRjGRADUATE’s performance does not degrade between SDR = R1 [ R2 does not holds necessarily true in all cases. the Kaggle DR detection test set, the most similar image- For instance, when an image is R3 by the international scale and grading-wise to the training data, and other datasets, due to the presence of > 20 HEMs, it would fit in the SDR cat- showing the generalization capability of the model. Further- egory of the Davis scale. However, through the conversion, this more, considering the misclassified cases, the 2 higher over- image would fall in the PPDR category, which can justify the diagnosis percentage of DRjGRADUATE in comparison to model’s confusion between R1 and R2 that is verified in Fig. 7e. the expected human performance moderately promotes patient follow-up, challenging specialists to reassess their opinion. De- Classification errors for R2 and R3 grades are mostly due to spite the high performance achieved, DRjGRADUATE tends the grading rules. Indeed, in cases where there are no other to fail to distinguish between R3 and R4 grades. As shown signs of R3 (Table 1), grading depends on HEM counting. in Figs. 15a and 15b, DRjGRADUATE learned to associate However, this task is non-trivial for specialists since HEM may photo-coagulation treatment and laser marks to R4, following appear in clusters, hindering lesion counting and thus introduc- the training set labels. Indeed, these images may not have signs ing noise on the data. Furthermore, the MIL assumption fails of R4, such as neovessels and pre-retinal hemorrhages, and thus for these cases, as there is no novel lesion to find, leading to pre- should either be marked as photocoagulated and removed from diction errors (Fig. 17b). Instead, DRjGRADUATE may have the dataset or be included with the grade corresponding to the associated lesion size to the grade. Alternatively, despite rely- visible lesions. This issue is further exacerbated by the high ing on a single lesion to perform grade inference, the construc- class imbalance of the data (Fig. 5), since the data augmentation tion of the lesion map L depends on the integration of features and class balancing scheme can contribute to the regulation of with di erent receptive fields and thus the lesion with the max- the weights update but not the limited variability of the samples. imum response may indirectly encode lesion counting, serving Ground truth R4 R3 R2 R1 R0 Ground truth R4 R3 R2 R1 R0 Ground truth R4 R3 R2 R1 R0 Ground truth R4 R3 R2 R1 R0 Ground truth R3 R2 R1 R0 Teresa Araujo ´ et al. / Medical Image Analysis (2020) 11 % images and the prediction. This allows the use of the common cat- 4 7 10 27 44 63 83 89 96 97 99 DMR egorical crossentropy as the main term of the loss function SCREEN-DR 5 9 14 28 42 57 73 84 95 97 99 and thus avoids the need for the direct optimization of the , Messidor-2 16 29 43 62 81 89 97 98 99 which requires higher computational power. On the other hand, 6 14 23 41 60 73 86 90 95 96 98 99 Kaggle-test DRjGRADUATE’s performance is inferior to the 0.84  re- 0.925 ported by Krause et al. (2018) on EyePACS-2 dataset (1818 images) (data publicly unavailable at the time of writing). How- 0.900 ever, their model was trained on an expensive dataset composed of over 1.6 million images, and with ground truth labels ad- 0.875 0.875 judicated by multiple ophthalmologists. Further, the authors 0.850 apply an empirically defined cascade of thresholds over their continuous network output in order to define the grade, ensur- 0.825 ing that a baseline sensitivity was achieved for the higher grade classes. Consequently, the system also tends to overgrade im- 0.800 ages, with a normalized over-grading of 72%, which can lead to a high number of unnecessary follow-ups. On the other hand, 0.775 DRjGRADUATE presents a lower over-estimation and thus a 0.750 more moderate trade-o between over and under diagnosis of DR. 0.725 5.2. Uncertainty estimation 0.700 0.15 0.21 0.29 0.4 0.54 0.7 0.89 1.14 Ideally, DRjGRADUATE should infer higher uncertainties uncertainty for misclassified cases. Indeed, high  values occur for low un- (a) Quadratic weighted Kappa vs uncertainty (xx-axis in logarithmic scale) for certainty (Fig. 8a), indicating that the DRjGRADUATE’s uncer- four di erent datasets (DMR: 4 classes). tainty estimation is a valid measure of the quality of the predic- Predicted 1.0 R0 R1 R2 R3 R4 tion. Fig. 8c further indicates that low uncertainty is associated 0.8 with correct or close classification, since the uncertainties are 0.27 0.35 0.39 0.43 0.38 considerably lower in the diagonal, and tend to increase with the 0.6 0.27 0.26 0.25 0.33 0.35 class distance. Furthermore, all grades follow a similar cumula- R0 0.4 0.31 0.31 0.25 0.25 0.28 tive histogram, as shown in Fig. 8b, indicating that the tendency R1 R2 illustrated in Fig. 8a does not result from the di erent propor- 0.38 0.41 0.27 0.20 0.23 0.2 R3 R4 tions of each grade among the analysed image subsets. Low 0.41 0.57 0.38 0.28 0.21 0.0 uncertainty of misclassifications tends to be due to the presence 0.5 1.0 1.5 uncertainty of lesions that greatly resemble another type, characteristic of (b) Cumulative percentage of images from (c) Average predicted uncer- a di erent grade (e.g., large MAs which are similar to HEMs), each grade included per uncertainty threshold, tainty matrix for the Kaggle or to the erroneous identification of acquisition errors as lesions for the Kaggle DR detection test set. DR detection test set. (Fig. 18). The lower  for higher uncertainty cases suggests that Fig. 8: Uncertainty of the predictions in relation with the DRjGRADUATE’s performance. DRjGRADUATE can be used cooperatively in the clinical prac- tice by asking ophthalmologists to review cases that are most likely wrong. Indeed, case selection via uncertainty may sig- as a proxy of the overall relevant lesion density (Fig. 17a). nificantly reduce the workload during the screening pipeline. For instance, discarding the 25% most uncertain images from 5.1.1. Comparison with the state-of-the-art Kaggle DR detection test set allows to achieve a 0.8 , a per- Overall, DRjGRADUATE achieves similar performance to formance similar to human experts (Krause et al., 2018), thus other state-of-the-art methods and shows high generalization allowing DRjGRADUATE to be used as an independent reader capability. Indeed, our model has a  similar to Takahashi capable of identifying dubious cases. et al. (2017)’s method on the DMR dataset (value computed High uncertainty tends to occur in very obscured and blurred by data provided by the authors) while being trained on a dif- images (Fig. 9a, 9b, 9c, 9e, 9f). In these cases, it is not possible ferent dataset and considering a di erent grading scale (inter- to make a diagnosis - at least not fully -, and thus predictions are national instead of Davis). Note that class-wise performance associated with high uncertainty estimations. This will be fur- comparison is limited because confusion matrices are not com- ther discussed in section 5.2.1. Fig. 9d shows an interesting case monly available in the literature. The proposed system achieves of a R4 image correctly predicted by DRjGRADUATE but with a performance similar to de la Torre et al. (2017), while having high uncertainty (0.999). This image shows a very large PFIB only 2% of the parameters and not requiring highly dedicated (PDR sign), occupying almost half of the FOV. However, as re- loss functions. Indeed, DRjGRADUATE’s architecture already ferred, DRjGRADUATE did not correctly learn these R4 signs promotes the distance minimization between the ground-truth (even less of this size) and thus it is not expected to be certain % images/class Quadratic weighted kappa Ground truth R4 R3 R2 R1 R0 12 Teresa Araujo ´ et al. / Medical Image Analysis (2020) (a) y ˆ = 1; y = 2. (b) y ˆ = 0; y = 0. (c) y ˆ = 2; y = 3. (d) y ˆ = 4; y = 4. (e) y ˆ = 2; y = 4. (f) y ˆ = 0; y = 2. g g g g g g (g) y ˆ = 4; y = 0. (h) y ˆ = 2; y = 2. (i) y ˆ = 3; y = 0. (j) y ˆ = 3; y = 4. (k) y ˆ = 1; y = 1. (l) y ˆ = 3; y = 3. g g g g g g (m) y ˆ = 2; y = 2. (n) y ˆ = 0; y = 0. (o) y ˆ = 2; y = 2. (p) y ˆ = 3; y = 2. (q) y ˆ = 3; y = 3. (r) y ˆ = 4; y = 4. g g g g g g Fig. 9: Uncertainty-associated predictions for images from the Kaggle DR detection test set. The color of the frame represents the uncertainty, y ˆ is the prediction and y the ground truth grade. Images were cropped around the field-of-view and resized to a square image. Uncertainty colorbar: 0:05 1:406 in making a diagnosis based on such image. This large white uncertainty (Fig. 9f). structure may have been interpreted as partially obscuring the However, since image quality is not trivial to define within image and thus the high uncertainty. Further, high uncertainties the context of clinical diagnosis, image quality assessment can also result from the resemblance of certain structures in the tends to be a subjective process (Wanderley et al., 2019). One images with DR signs, such as CWSs (Fig. 9i), and with pho- can consider more than two levels of quality, such as good, par- tocoagulation marks (Fig. 9g). Contrastingly, DRjGRADUATE tial, bad, where a partial quality image still has a region suitable correctly predicts, with low uncertainty, the good quality and for diagnosis. However, most datasets are labeled in a binary diagnosable images from Fig. 9k-9r, apart from Fig. 9p. We be- manner (good or bad quality), forcing partially obscured im- lieve this case is an example of the noise present in Kaggle DR ages to considered bad or good. detection dataset annotations, since the image is labeled as R2, Further, as image quality is not the only factor influencing but shows several CWSs (R3 sign), leading DRjGRADUATE to the uncertainty of the model prediction, a full correspondence predict R3 with very low uncertainty (0.05). between image quality and the uncertainty values is not likely. However, a tendency of bad quality images being associated, overall, with higher uncertainties is expected. 5.2.1. Low quality images Poor image quality can lead to the occlusion or distortion DRjGRADUATE produced significantly higher uncertainties of lesions, which hinders a meaningful and complete clinical for bad than for good quality images in each dataset (Fig. 10a), decision (Niemeijer et al., 2006). Lack of focus, sharpness or which suggests it can potentially be used in the DR diagnosis illumination often characterize bad quality images (Pires Dias pipeline without the need for an a priori image quality detec- et al., 2014). In theory, and as suggested by Fig. 9, the algorithm tion step. The HRF datasets is particularly interesting since it should output higher uncertainties for images of poor quality, contains pairs of good/bad quality images, but has the incon- since in the these cases it can not be as sure of the diagnosis venient of having only 18 pairs of images. Some of the HRF as in a good quality image. Indeed, Fig. 9a-f mainly show ex- bad quality images di er only slightly from the good quality amples of largely obscured images which are associated with counterparts (Fig. 19), which leads to similar uncertainty esti- very high uncertainty estimations. Usually, when the image is mation between the two sets. Overall, the bad quality images partially available for diagnosis, if DRjGRADUATE can find seem to still be suitable for diagnosis, which may justify the lesions in the viable image region it predicts a disease grade lower uncertainties for the bad quality subset comparing with with high uncertainty (HEMs from R2, in Fig. 9h), whereas if what is verified on the other datasets. DRIMDB is the more it cannot find any signs it predicts it as healthy, again with high dichotomic dataset, where images are clearly bad or good, ap- Teresa Araujo ´ et al. / Medical Image Analysis (2020) 13 p < 0.05 p < 0.05 d ≈ 0.90 d ≈ 2.79 p < 0.05 d ≈ 0.68 p < 0.05 d ≈ 0.77 p < 0.05 1.2 p < 0.05 d ≈ 0.66 d ≈ 1.48 0.8 p < 0.05 d ≈ 0.85 p < 0.05 0.8 p < 0.05 p > 0.05 p > 0.05 d ≈ 0.09 d ≈ 0.51 d ≈ 0.18 d ≈ 0.05 0.4 0.4 0.0 0.0 Kaggle test DR1 Giana ISIC Cataracts BACH Screen-DR dmr IDRID Messidor-2 HRF DRIMDB Datasets Datasets (a) Uncertainties predicted by DRjGRADUATE for good and bad (b) Uncertainties predicted by DRjGRADUATE for unfamiliar data types from di erent medical quality images from public datasets. p-values obtained through image datasets and for the studied eye fundus images datasets. p-values obtained through the the Kruskal-Wallis H Test and Cohen’s d between each good and Kruskal-Wallis H Test and Cohen’s d between each dataset and its right neighbour on the plot bad sets from a given dataset are presented. are presented. Fig. 10: Uncertainty prediction for di erent types of tasks. mentation scheme performed during DRjGRADUATE’s train- ing, which contemplated brightness adjustments and contrast normalization, and most likely made the network robust to these changes. As the major factor accounting for bad quality in the tests with the bad/good quality datasets (Section 3.1.2 (a) std = 3 (b) std = 6 (c) std = 10 (d) std = 20 and Fig. 10a) was blur/decreased sharpness, not image con- trast/illumination conditions, the sensitivity analysis results are Fig. 11: Examples of convolution with Gaussian filters of increasing standard according to the expectations. deviation (std). 5.2.2. Unfamiliar data types The uncertainty estimation on the unfamiliar data types sug- 1.75 gests that unfamiliar medical image types, i.e., which are not 1.50 eye fundus images, could be detected as outliers by the analysis 1.25 of the uncertainty values produced by the algorithm, which are 1.00 in average higher than for the eye fundus images. 0.75 Colonoscopy images are probably the most similar to the eye 0.50 fundus, showing a similar circular FOV and color (Fig. 10b), 0.25 and sometimes even red and bright structures which may re- 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 semble retinal lesions and anatomical structures. Considering std that the DRjGRADUATE was trained with eye fundus images, Fig. 12: Uncertainty estimation of DRjGRADUATE for images resulting from images that most resemble these can lead to more uncertainty convolution with Gaussian filters of increasing standard deviation (std). The since they are somehow familiar but do not present the usual uncertainty values are the average of the uncertainties for 5 datasets (Kaggle and expected structures for diagnosis, similarly to what happens DR detection test set, Messidor-2, IDRID, SCREEN-DR and DMR) consider- ing a 95% confidence interval. in bad quality eye fundus images. H&E stained microscopy im- ages of breast tissue are the closest to the eye fundus in terms of uncertainty estimation, and are visually far from retinal images. This may be due to the presence of several small rounded struc- parently not presenting a middle term. This is reflected on the tures, darker than the background (cells’ nuclei) which may estimated uncertainties, with the bad quality set showing an av- resemble retinal lesions. Considering that our MIL approach erage uncertainty of 0.750.22 and the good set of 0.300.12. learns to identify the local structures responsible for a given DR Note that both HRF and DR1 contemplate only blur in the bad grade, it makes sense that resemblance at a local level (as hap- quality images, and thus low contrast or bad illumination are pens in microscopy images) leads to lower uncertainties than not considered. resemblance at the global image level (colonoscopy images). Regarding the sensitivity analysis with the blurring opera- tions, we can see based on Fig. 12 that the uncertainty shows a 5.3. Explainability tendency to increase with the increase of the image blur, which decreases image quality. This corroborates our previous find- DRjGRADUATE’s MIL-based architecture leads to the ex- ings that image quality is correlated with uncertainty. When it planation maps to pinpoint regions relevant for diagnosis, and comes to the contrast variation experiments, we believe that the thus highlighting every single lesion is not expected. Because uncertainty did not show any correlation with it due to data aug- of this, the global lesion-wise is not reliable enough to fully Uncertainty uncertainty Uncertainty 14 Teresa Araujo ´ et al. / Medical Image Analysis (2020) Explanation map Annotated lesions Explanation map Annotated lesions MA MA (a) y ˆ = 1; y = 1. (b) y ˆ = 0; y = 1. HEM HEM HEM MA MA (c) y ˆ = 2; y = 2. (d) y ˆ = 3; y = 2. g g CWS HEM IrMA HEM HEM (f) y ˆ = 3; y = 3. (e) y ˆ = 3; y = 3. g NV PHEM NV (g) y ˆ = 4; y = 4. (h) y ˆ = 3; y = 4. Fig. 13: Maps predicted by our method along with the ophthalmologist’s annotated lesions in the SCREEN-DR dataset. The curves are the contours of the explanation maps (threshold=0.3), and the red square indicates the region of most relevance for diagnosis (corresponding to the maximum in the network’s output activation map). Close-ups of relevant regions are shown, with 3 indicating that the region is correctly predicted, and 7 indicating otherwise. Below each close-up the name of the corresponding ground truth lesion (if existent) is shown. Explanation map:  R1,  R2,  R3,  R4; ground truth map:  MA,  HEM,  EX, CWS, IrMA,  NV,  PHEM, PFIB. characterize an image. The detection performance is further that the model detects this lesion type both as a R2 or a R3 hindered by the fact that specialists may not annotate all ab- sign (Fig. 13d), which is understandable since a low number of normalities (Fig. 13) since this can be a large time consuming HEMs is a sign for R2 and a high number a sign of R3. Despite task. On the other hand, DRjGRADUATE is able to correctly this image being annotated as R2, the high number of HEMs identify plausible regions for the reasoning behind the grade of may justify the R3 prediction. However, DRjGRADUATE fails a given image. The analysis of the produced attention maps to detect NVs and other R4 lesions such as PHEMs (Fig. 13g) show the ability of the system to detect lesions as small as MAs and PFIB (Fig. 13h). Most of the times the reason for a R4 pre- (Fig. 13a). As expected, when image quality is poor, image diction is related with the presence of photocoagulation marks diagnosis is hindered and lesions are often missed (Fig. 13b). (Fig. 13g) which corroborates the idea that these were learn by Confusion between MAs and HEMs occurs for some images the algorithm as being R4 signs. Nevertheless, R4 lesions are (Fig. 13c), which is justifiable by the high resemblance of small sometimes highlighted in the maps as signs of disease. For in- HEMs and MAs. The detection of R3 is particularly interest- stance, NVs and PHEMs are detected as R2 signs in these two ing due to the variety of configurations that are associated with images, since they slightly resemble HEMs. this grade (Table 1), which DRjGRADUATE seems to be able Interestingly, the inherent explainability of DRjGRADUATE to identify: a high number of HEMs (Fig. 13d), some CWSs allows for the a posteriori identification of systematic errors (Fig. 13e) or IrMAs (Fig. 13f). Regarding HEMs, one can see during image acquisition. For instance, several Messidor-2 im- Teresa Araujo ´ et al. / Medical Image Analysis (2020) 15 0.80 Omax 0.75 Oclass 0.70 0.65 0.60 0.55 0.50 0.45 0.1 0.2 0.3 0.4 0.5 0.6 0.7 threshold Fig. 14: O and O for di erent threshold values (0.1, 0.3, 0.5 and 0.7) max class (b) y ˆ = 4; y = 4 (Kaggle DR detec- (a) y ˆ = 4; y = 4 (Messidor-2). applied to the attention maps (Eq. 3). tion test). ages have been wrongly classified as R1 due to the presence of dust particles in the camera lens (Fig.18). However, the clear explanation provided would easily allow the clinician to correct the acquisition error and thus improve the overall quality of the screening pipeline. 6. Conclusions (d) y ˆ = 4; y = 0 (Kaggle DR detec- (c) y ˆ = 4; y = 0 (Messidor-2). tion test). Fig. 15: Examples of images from Messidor-2 and Kaggle DR detection test We propose DRjGRADUATE, a novel deep learning-based sets that are predicted by DRjGRADUATE as R4. The curves are the contours approach for DR grading that provides an uncertainty and ex- of the attention maps (threshold=0.3).  R1, R2, R3,  R4. planation associated with each prediction. DRjGRADUATE was trained on the Kaggle DR detection training set, and Acknowledgments achieved state-of-the-art performance in several DR-labeled datasets. The classification of the R4 grade was the least accom- Teresa Araujo ´ is funded by the FCT grant contract plished task, suggesting the training did not correctly capture SFRH/BD/122365/2016. Guilherme Aresta is funded by the the PDR signs. This is most likely due to the type of annotation FCT grant contract SFRH/BD/120435/2016. This work is fi- in the Kaggle DR detection dataset, in which images that under- nanced by the ERDF European Regional Development Fund went photocoagulation treatment and show laser marks appear through the Operational Programme for Competitiveness and to be labeled as PDR even when not presenting the typical signs Internationalisation - COMPETE 2020 Programme, and by Na- of this grade. tional Funds through the FCT - Fundac ¸ao ˜ para a Ciencia ˆ e a The analysis of the predictions’ uncertainties showed the ten- Tecnologia within project CMUP-ERI/TIC/0028/2014. dency of lower  for higher uncertainty cases, thus suggesting that the estimated uncertainty is a viable measure of the quality of the prediction. This means that at test time, i.e. where the References performance of the system is not known, the uncertainty can be Abramo ` , M.D., Lou, Y., Erginay, A., Clarida, W., Amelon, R., Folk, J.C., a measure of how much a specialist can trust the algorithm’s Niemeijer, M., 2016. Improved Automated Detection of Diabetic Retinopa- decision, allowing to rank the predictions and select which im- thy on a Publicly Available Dataset Through Integration of Deep Learn- ages would benefit most from a second evaluation. Further, the ing. Investigative Opthalmology & Visual Science 57, 5200. doi:10.1167/ explanation map produced by our algorithm should be easily iovs.16-19964. interpretable by the retinal specialists, and thus constitutes a Abramo , M.D., Suttorp-Schulten, M.S., 2005. Web-based screening for di- abetic retinopathy in a primary care population: The EyeCheck Project. valuable tool for mitigating the black-box behaviour associated Telemedicine Journal and e-Health 11, 668–674. doi:10.1089/tmj.2005. with deep NNs. 11.668. Acharya, U.R., Ng, E.Y.K., Tan, J.H., Sree, S.V., Ng, K.H., 2012. An inte- Comparing with the state-of-the-art, DRjGRADUATE has grated index for the identification of diabetic retinopathy stages using tex- the advantage of producing both an uncertainty estimation asso- ture parameters. Journal of Medical Systems 36, 2011–2020. doi:10.1007/ ciated with the prediction and an explanation map highlighting s10916-011-9663-8. the most relevant regions for the classification, without compro- Arev ´ alo, J.F., Lasave, A.F., Zeballos, D.G., Bonafonte-Royo, S., 2013. Diabetic Retinopathy, in: Retinal and Choroidal Manifestations of Selected Systemic mising the DR grading performance. We believe that coupling Diseases. Springer New York, New York, NY, pp. 387–416. doi:10.1007/ an uncertainty estimation and explanation with the DR grade 978-1-4614-3646-1{\_}21. will ease the application of DRjGRADUATE in a clinical set- Bach, S., Binder, A., Montavon, G., Klauschen, F., Muller ¨ , K.R., Samek, W., ting. 2015. On pixel-wise explanations for non-linear classifier decisions by layer- hits 16 Teresa Araujo ´ et al. / Medical Image Analysis (2020) (a) y ˆ = 0; y = 1. (b) y ˆ = 1; y = 0. (a) y ˆ = 3; y = 3. (b) y ˆ = 3; y = 2. g g g g Fig. 17: Examples of R2 and R3 Kaggle DR detection dataset images (not used for training) along with DRjGRADUATE predictions. The curves are the contours of the attention maps (threshold=0.3).  R1, R2, R3, R4. (c) y ˆ = 1; y = 2. (d) y ˆ = 1; y = 2. g g Fig. 16: Examples of R0-R2 Kaggle DR detection dataset images (not used for training) along with the DRjGRADUATE predictions. The curves are the contours of the attention maps (threshold=0.3).  R1, R2, R3, R4. (a) y ˆ = 1; y = 0. (b) y ˆ = 1; y = 0. g g Fig. 18: Examples of healthy Messidor-2 images mispredicted as R1. These wise relevance propagation. PLoS ONE 10, 1–46. doi:10.1371/journal. images show camera artifacts (present in the same location across several dif- pone.0130140. ferent images) which are detected as R1 typical lesions. The curves are the Buda, M., Maki, A., Mazurowski, M.A., 2018. A systematic study of the class contours of the attention maps (threshold=0.3).  R1,  R2, R3, R4. imbalance problem in convolutional neural networks. Neural Networks 106, 249–259. doi:10.1016/j.neunet.2018.07.011. Costa, P., Araujo, ´ T., Aresta, G., Galdran, A., Mendonc ¸a, A.M., Smailagic, A., Campilho, A., 2019. EyeWeS : Weakly Supervised Pre-Trained Convolu- Gupta, G., Kulasekaran, S., Ram, K., Joshi, N., Sivaprakasam, M., Gandhi, R., tional Neural Networks for Diabetic Retinopathy Detection, in: International 2017. Local characterization of neovascularization and identification of pro- Conference on Machine Vision Applications (MVA), Tokyo. pp. 5–15. liferative diabetic retinopathy in retinal fundus images. Computerized Med- Cuadros, J., Bresnick, G., 2009. EyePACS: An adaptable telemedicine system ical Imaging and Graphics 55, 124–132. doi:10.1016/j.compmedimag. for diabetic retinopathy screening. Journal of Diabetes Science and Tech- 2016.08.005. nology 3, 509–516. doi:10.1177/193229680900300315. Hall, A.B., 2011. Recognising and managing diabetic retinopathy How I look Dutra Medeiros, M., Mesquita, E., Papoila, A.L., Genro, V., Raposo, J.F., 2015. for diabetic retinopathy : a vision technician s experience. Community Eye First diabetic retinopathy prevalence study in Portugal: RETINODIAB Health 24. StudyEvaluation of the screening programme for Lisbon and Tagus Valley Harangi, B., Hajdu, A., 2013. Improving automatic exudate detection based region. British Journal of Ophthalmology 99, 1328–1333. doi:10.1136/ on the fusion of the results of multiple active contours. Proceedings - Inter- bjophthalmol-2015-306727. national Symposium on Biomedical Imaging , 45–48doi:10.1109/ISBI. Friston, K., 2010. The free-energy principle: A unified brain theory? Nature 2013.6556408. Reviews Neuroscience 11, 127–138. doi:10.1038/nrn2787. Kendall, A., Gal, Y., 2017. What uncertainties do we need in Bayesian deep Gal, Y., Ghahramani, Z., 2016. Dropout as a Bayesian Approximation: Rep- learning for computer vision?, in: Advances in Neural Information Process- resenting Model Uncertainty in Deep Learning, in: Proceedings of the 33rd ing Systems, Long Beach. pp. 5575–5585. International Conference on Machine Learning, New York. doi:10.1109/ Kohler ¨ , T., Budai, A., Kraus, M.F., Odstrcilik, ˇ J., Michelson, G., Hornegger, J., TKDE.2015.2507132. 2013. Automatic no-reference quality assessment for retinal fundus images Garg, S., Davis, R.M., 2009. Diabetic retinopathy screening update. Clinical using vessel segmentation, in: Proceedings of CBMS 2013 - 26th IEEE In- Diabetes 27, 140–145. doi:10.2337/diaclin.27.4.140. ternational Symposium on Computer-Based Medical Systems, pp. 95–100. Gargeya, R., Leng, T., 2017. Automated Identification of Diabetic Retinopa- doi:10.1109/CBMS.2013.6627771. thy Using Deep Learning. Ophthalmology , 1–8doi:10.1016/j.ophtha. Krause, J., Gulshan, V., Rahimy, E., Karth, P., Widner, K., Corrado, G.S., Peng, 2017.02.008. L., Webster, D.R., 2018. Grader Variability and the Importance of Reference Gonzalez-Gonzalo, ´ C., Liefers, B., Ginneken, B.V., Sanchez, ´ C.I., 2018. Im- Standards for Evaluating Machine Learning Models for Diabetic Retinopa- proving weakly-supervised lesion localization with iterative saliency map thy. Ophthalmology 125, 1264–1272. doi:10.1016/j.ophtha.2018.01. refinement , 4–6. 034. Gulshan, V., Peng, L., Coram, M., Stumpe, M.C., Wu, D., Narayanaswamy, Kruskal, W.H., Wallis, W.A., Kruskal, W.H., Wallis, W.A., 2012. Journal of the A., Venugopalan, S., Widner, K., Madams, T., Cuadros, J., Kim, R., American Statistical Analysis. Journal of the American Statistical Associa- Raman, R., Nelson, P.C., Mega, J.L., Webster, D.R., 2016. Develop- tion 47, 583–621. doi:10.2307/2280779. ment and Validation of a Deep Learning Algorithm for Detection of Di- Leibig, C., Allken, V., Ayhan, M.S., Berens, P., Wahl, S., 2017. Leveraging abetic Retinopathy in Retinal Fundus Photographs. Jama 304, 649–656. uncertainty information from deep neural networks for disease detection. doi:10.1001/jama.2016.17216. Scientific Reports 7, 1–14. doi:10.1038/s41598-017-17876-z. Teresa Araujo ´ et al. / Medical Image Analysis (2020) 17 jgme-d-12-00156.1. Takahashi, H., Tampo, H., Arai, Y., Inoue, Y., Kawashima, H., 2017. Applying artificial intelligence to disease staging: Deep learning for improved staging of diabetic retinopathy. Plos One 12, e0179790–e0179790. doi:10.1371/ journal.pone.0179790. The Royal College of Ophthalmologists, 2012. Diabetic Retinopathy Guide- lines. Technical Report December. de la Torre, J., Puig, D., Valls, A., 2018. Weighted kappa loss function for multi- class classification of ordinal data in deep learning. Pattern Recognition Letters 105, 144–154. doi:10.1016/j.patrec.2017.05.018. de la Torre, J., Valls, A., Puig, D., 2017. A Deep Learning Interpretable Classi- fier for Diabetic Retinopathy Disease Grading. arXiv . Wanderley, D.S., Araujo, T., Carvalho, C.B., Maia, C., Penas, S., Carneiro, A., (a) Bad quality image, u=0.415. (b) Good quality image, u=0.295. Mendonca, A.M., Campilho, A., 2019. Analysis of the performance of spe- cialists and an automatic algorithm in retinal image quality assessment, in: 2019 IEEE 6th Portuguese Meeting on Bioengineering (ENBENG), IEEE. Fig. 19: Example of a pair of good/bad quality images from the HRF dataset. doi:10.1109/ENBENG.2019.8692506. Wild, 2004. Estimates for the year 2000 and projections for 2030. World Health 27, 1047–1053. doi:10.2337/diacare.27.5. Mahendran, G., Dhanasekaran, R., 2015. Investigation of the severity level 1047DiabetesCareMay2004vol.27no.51047-1053. of diabetic retinopathy using supervised classifier algorithms. Computers Wu, L., Fernandez-Loaiza, P., Sauma, J., Hernandez-Bogantes, E., Masis, M., and Electrical Engineering 45, 312–323. doi:10.1016/j.compeleceng. 2013. Classification of diabetic retinopathy and diabetic macular edema. 2015.01.013. World journal of diabetes 4, 290–4. doi:10.4239/wjd.v4.i6.290. Massin, P., Chabouis, A., Erginay, A., Viens-Bitker, C., Lecleire-Collet, A., Zeiler, M.D., Fergus, R., 2014. Visualizing and Understanding Convolutional Meas, T., Guillausseau, P.J., Choupot, G., Andre, ´ B., Denormandie, P., Networks, in: Fleet, D., Pajdla, T., Schiele, B., Tuytelaars, T. (Eds.), Com- 2008. OPHDIAT c : A telemedical network screening system for diabetic puter Vision – ECCV 2014, Springer International Publishing, Cham. pp. retinopathy in the Ile-de-France. Diabetes and Metabolism 34, 227–234. 818–833. doi:10.1016/j.diabet.2007.12.006. Zhang, B., Wu, X., You, J., Li, Q., Karray, F., 2010. Detection of microa- Nguyen, T.T., Wong, T.Y., 2009. Retinal vascular changes and diabetic neurysms using multi-scale correlation coecients. Pattern Recognition 43, retinopathy. Current Diabetes Reports 9, 277–283. doi:10.1007/ 2237–2248. s11892-009-0043-4. Niemeijer, M., Abramo ` , M.D., van Ginneken, B., 2006. Image structure clustering for image quality verification of color retina images in diabetic retinopathy screening. Medical Image Analysis 10, 888–898. doi:10.1016/ j.media.2006.09.006. Pires, R., Jelinek, H.F., Wainer, J., Rocha, A., 2012. Retinal image quality analysis for automatic diabetic retinopathy detection. Brazilian Sympo- sium of Computer Graphic and Image Processing , 229–236doi:10.1109/ SIBGRAPI.2012.39. Pires Dias, J.M., Oliveira, C.M., Da Silva Cruz, L.A., 2014. Retinal image qual- ity assessment using generic image quality indicators. Information Fusion 19, 73–90. doi:10.1016/j.inffus.2012.08.001. Porwal, P., Pachade, S., Kamble, R., Kokare, M., Deshmukh, G., Sahasrabud- dhe, V., Meriaudeau, F., 2018. Indian Diabetic Retinopathy Image Dataset (IDRiD). IEEE Dataport . Prentasic, P., Loncaric, S., 2016. Detection of exudates in fundus pho- tographs using deep neural networks and anatomical landmark detection fusion. computer methods and programs in biomedicine 137, 281–292. doi:10.1109/ISPA.2015.7306056. Quellec, G., Charriere, ` K., Boudi, Y., Cochener, B., Lamard, M., 2017. Deep image mining for diabetic retinopathy screening. Medical Image Analysis 39, 178–193. doi:10.1016/j.media.2017.04.012. Sawilowsky, S.S., 2009. New E ect Size Rules of Thumb. Journal of Modern Applied Statistical Methods 8, 597–599. doi:10.22237/jmasm/ Schully, T., 2012. Diabetes in Numbers. Nature 485, S2–S3. doi:10.1038/ nn.2369. Scotland, G.S., McNamee, P., Philip, S., Fleming, A.D., Goatman, K.A., Prescott, G.J., Fonseca, S., Sharp, P.F., Olson, J.A., 2007. Cost-e ectiveness of implementing automated grading within the national screening pro- gramme for diabetic retinopathy in Scotland. The British journal of oph- thalmology 91, 1518–23. doi:10.1136/bjo.2007.120972. Sevik, U., Kose, ¨ C., Berber, T., Erdol, ¨ H., 2014. Identification of suitable fundus images using automated quality assessment methods. Journal of Biomedical Optics 19, 046006. doi:10.1117/1.jbo.19.4.046006. Sil Kar, S., Maity, S., 2017. Automatic Detection of Retinal Lesions for Screen- ing of Diabetic Retinopathy. IEEE Transactions on Biomedical Engineering 9294, 1–1. doi:10.1109/TBME.2017.2707578. Simonyan, K., Vedaldi, A., Zisserman, A., 2014. Deep Inside Convolutional Networks: Visualising Image Classification Models and Saliency Maps, in: International Conference on Learning Representations (ICRL), Calgary. Sullivan, G.M., Feinn, R., 2012. Using E ect Sizeor Why the P Value Is Not Enough . Journal of Graduate Medical Education 4, 279–282. doi:10.4300/ http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Electrical Engineering and Systems Science arXiv (Cornell University)

DR$\vert$GRADUATE: uncertainty-aware deep learning-based diabetic retinopathy grading in eye fundus images

Loading next page...
 
/lp/arxiv-cornell-university/dr-vert-graduate-uncertainty-aware-deep-learning-based-diabetic-dZMWRCMBNx

References (55)

ISSN
1361-8415
eISSN
ARCH-3348
DOI
10.1016/j.media.2020.101715
Publisher site
See Article on Publisher Site

Abstract

Medical Image Analysis (2020) Contents lists available at ScienceDirect Medical Image Analysis journal homepage: www.elsevier.com/locate/media DRjGRADUATE: uncertainty-aware deep learning-based diabetic retinopathy grading in eye fundus images a,b a,b c d,e d d,e Teresa Araujo ´ , Guilherme Aresta , Lu´ ıs Mendonc ¸a , Susana Penas , Carolina Maia , Angela Carneiro , Ana Maria a,b a,b Mendonc ¸a , Aurelio ´ Campilho INESC TEC - Institute for Systems and Computer Engineering, Technology and Science, 4200-465 Porto, Portugal Faculty of Engineering of University of Porto, 4200-465 Porto, Portugal Department of Ophthalmology of Hospital de Braga, Braga, Portugal Department of Ophthalmology of Centro Hospitalar Sao ˜ Joao, ˜ Porto, Portugal Department of Surgery and Physiology of Faculty of Medicine of University of Porto, Porto, Portugal A R T I C L E I N F O A B S T R A C T Diabetic retinopathy (DR) grading is crucial in determining the adequate treatment Article history: and follow up of patients, but the screening process can be tiresome and prone to er- rors. Deep learning approaches have shown promising performance as computer-aided diagnosis (CAD) systems, but their black-box behaviour hinders clinical application. Diabetic retinopathy grading, Deep learning, Uncertainty, Explainability We propose DRjGRADUATE, a novel deep learning-based DR grading CAD system that supports its decision by providing a medically interpretable explanation and an es- timation of how uncertain that prediction is, allowing the ophthalmologist to measure how much that decision should be trusted. We designed DRjGRADUATE taking into account the ordinal nature of the DR grading problem. A novel Gaussian-sampling ap- proach built upon a Multiple Instance Learning framework allow DRjGRADUATE to infer an image grade associated with an explanation map and a prediction uncertainty while being trained only with image-wise labels. DRjGRADUATE was trained on the Kaggle DR detection training set and evaluated across multiple datasets. In DR grading, a quadratic-weighted Cohen’s kappa () between 0.71 and 0.84 was achieved in five dif- ferent datasets. We show that high  values occur for images with low prediction uncer- tainty, thus indicating that this uncertainty is a valid measure of the predictions’ quality. Further, bad quality images are generally associated with higher uncertainties, showing that images not suitable for diagnosis indeed lead to less trustworthy predictions. Ad- ditionally, tests on unfamiliar medical image data types suggest that DRjGRADUATE allows outlier detection. The attention maps generally highlight regions of interest for diagnosis. These results show the great potential of DRjGRADUATE as a second- opinion system in DR severity grading. c 2020 Elsevier B. V. All rights reserved. c 2019. Licensed under the Creative Commons CC-BY-NC-ND 4.0 li- cense http://creativecommons.org/licenses/by-nc-nd/4.0/. guilherme.m.aresta@inesctec.pt (Guilherme Aresta), Corresponding author amendon@fe.up.pt (Ana Maria Mendonc ¸a), campilho@fe.up.pt (Aurelio ´ e-mail: tfaraujo@inesctec.pt (Teresa Araujo), ´ Campilho) arXiv:1910.11777v2 [eess.IV] 29 May 2020 2 Teresa Araujo ´ et al. / Medical Image Analysis (2020) A B C Table 1: International dDiabetic retinopathy (DR) severity scale. NPDR: Non- proliferative DR; PDR: proliferative DR. Non-referable DR: R0 and R1; Refer- able DR: R2, R3 and R4. Microaneurysms Hemorrhages Hard exudates Grade Description R0 - No DR No visible sign of abnormalities Mild R1 - Presence of microaneurysms only NPDR Moderate More than just microaneurysms but Cotton-wool Neovascu- R2 - larization spots NPDR less than severe NPDR Any of the following: Fig. 1: Examples of typical retinal lesions characteristic of diabetic retinopathy. >20 intraretinal hemorrhages Severe R3 - Venous beading NPDR Intraretinal microvascular abnormalities no signs of PDR 1. Introduction Either or both of the following: R4 - PDR Neovascularization Diabetic retinopathy (DR) is a complication of diabetes and Vitreous/pre-retinal hemorrhage is one of the leading causes of blindness worldwide (The Royal College of Ophthalmologists, 2012; Wild, 2004) and the num- ber of diabetic patients is expected to grow from 346 million in 2012 to 552 million people in 2030 (Schully, 2012; The Royal DR is generally classified according to the absence or pres- College of Ophthalmologists, 2012). The majority of visual ence of new abnormal vessels as non-proliferative DR (NPDR) loss cases can be prevented with early detection and adequate or proliferative DR (PDR), respectively. NPDR can be further treatment (Garg and Davis, 2009; Wu et al., 2013). However, graded as mild, moderate and severe (The Royal College of one third of the diabetic patients su er from DR without hav- Ophthalmologists, 2012). These stages can be ordered accord- ing symptoms of visual problems, leading to the progression of ing to their risk of progression (Arev ´ alo et al., 2013). Namely, the disease without treatment (Acharya et al., 2012). Conse- Table 1 represents the International Clinical DR Scale (?), quently, regular check-ups via screening programs are essential adapted to consider only one image for diagnosis (?) by remov- to achieve early diagnosis of DR. ing the quadrant concept, since in the case of a single fundus image analysis quadrant assessment is limited due to the under- DR screening programs have been developed across 1 2 representation of some of the quadrants. Thus, in practice, the the world, namely in the United Kingdom , Ireland exact counting of lesions in each quadrant is not performed by Netherlands (Abramo and Suttorp-Schulten, 2005), 4 5 the specialists, but rather the estimation of what happens in the France (Massin et al., 2008), United States (Cuadros full eye based on their experience and medical knowledge. and Bresnick, 2009), Canada (Nguyen and Wong, 2009) and Portugal (Dutra Medeiros et al., 2015). In these programs, DR Mild NPDR is characterized by the presence of microa- detection and DR severity grading are performed by trained neurysms (MAs), which are the first sign of DR. Hemorrhages specialists by visual analysis of color eye fundus images. (HEMs) start to appear with the evolution of the disease in the These images are captured by fundus photography, allowing deeper layers of the retina. The disease progression results in an for non-invasive visual assessment of several retinal structures, increase of MAs, dot-hemorrhages (red lesions), retinal thick- as depicted in Fig. 1: 1) the fovea, a reddish area without ves- ening, hard exudates (EXs), and the appearance of intraretinal sels located at the center of the macula; 2) the optic disc (OD) microvascular abnormalities (IrMAs), venous beading and cot- which is a round structure close to the macula that corresponds ton-wool spots (CWSs), also known as soft exudates, in the in- to the exit of the axons of ganglion cells that constitute the nermost retina. Finally, PDR results from abnormal vessel for- optic nerve and the door for the retinal vessels that supply the mation, i.e., neovessels (NVs). Other signs of PDR are the vit- retina, and 3) blood vessels, that radiate from the OD center. reous/pre-retinal hemorrhages (PHEMs) and pre-retinal fibrosis During screening, specialists search for abnormalities in these (PFIB) (Hall, 2011). Fig. 1 shows examples of some of these structures and in the retinal tissue and classify the severity of lesions. the disease according to the type and quantity of the findings. Computer-aided diagnosis (CAD) systems can improve the Despite the overall success of these programs, the diagnosis DR screening pipeline both reducing the burden (Scotland et al., process is still prone to errors due to the large number of 2007) and by proving a second objective opinion to the oph- patients, poor image acquisition and variety of lesions. thalmologists, reducing subjectivity in the diagnosis (Gargeya and Leng, 2017; Abramo ` et al., 2016; Gulshan et al., 2016; Quellec et al., 2017). Namely, deep learning has recently al- https://www.gov.uk/topic/population-screening-programmes/ lowed for these systems to achieve near-human performance in diabetic-eye DR detection, i.e. detection of at least one DR-related lesion. https://www.diabeticretinascreen.ie On the other hand, DR grading, i.e. staging of the pathology http://www.eyecheckklant.nl according to its severity, is a more complex problem since it re- reseau-ophdiat.aphp.fr https://www.eyepacs.com quires the identification and integration of di erent abnormal- Teresa Araujo ´ et al. / Medical Image Analysis (2020) 3 ities. In fact, specialists tend to disagree on the grade of com- performed by a 11 convolution and global max-pooling, fol- plex cases (Krause et al., 2018) and thus the external opinion lowing the assumption that the presence of one unhealthy patch of CAD systems may contribute to further improve the success is enough to label the entire image. Because of this, the out- of DR grading. However, due to this higher complexity and the put of this convolution is a map where each element indicates inherent black box behaviour of deep learning systems, it may the malignancy score of a corresponding patch of the input im- be dicult for specialists to understand and accept the system’s age, thus explaining why the image label was predicted. In- decision. With this in mind, we propose a novel DR grading stead, Gonzalez-Gonzalo ´ et al. (2018) has applied an iterative CAD system that not only supports its decision by providing a inpainting approach that, at each step, allows to identify pro- medically interpretable explanation but also estimates how un- gressively less discriminative pixels of the image for DR classi- certain that prediction is. By doing so, the ophthalmologist not fication. This method involves the computation of the saliency only has a reasoning behind the system’s decision, but also a map at each step, i.e., the derivative of the classification output measure of how much that decision should be trusted, making with respect to the input image using guided back-propagation. the DR diagnosis process less prone to errors. The authors use a pre-trained VGG-16 on the Kaggle dataset (512512 pixel images). Similarly to other grading approaches, the problem is treated as being ordinal and thus the model is op- 1.1. State-of-the-art timized using a mean-square error loss. To avoid overfitting due Deep learning has become the preferred approach, over field to the high class imbalance, the classes are balanced taking into knowledge-based feature engineering approaches, for DR grad- account the referal/non-referal sample distribution. ing since it allows to better exploit the large number of data Besides explaining the decision, uncertainty estimation of available and to better deal with the labeling noise resulting a model’s prediction is of interest as it helps to avoid conse- from the complexity of the task. For instance, Krause et al. quences of the blind use of the model response (Kendall and (2018) fine-tuned an Inception-v4 with 1,600,000 images of Gal, 2017). This is particularly true on medical settings, where size 779779 pixels to avoid the elimination of small lesions misdiagnosis can have severe consequences for the patients. such as MAs. Also, to account for data imbalance, the authors With this in mind, Leibig et al. (2017) analysed the uncer- apply a cascade of thresholds on the output probabilities of the tainty information from deep neural nets for DR detection and model for each DR grade to obtain the final image classification. referable DR detection. The authors tested the dropout-based These thresholds are sequentially applied in a severity descend- Bayesian uncertainty estimation, comparing it with alternative ing order to ensure that a minimum grade-wise sensitivity is techniques, such as directly analysing the network’s softmax achieved. Instead, de la Torre et al. (2018) considered DR grad- output. In the dropout-based approach, dropout is switched ing as an ordinal classification problem and thus aimed at reduc- on during inference. Performing multiple predictions on the ing the distance between predicted and true labels of an image. same input image can be interpreted as having an ensemble of For that, they proposed a quadratic-weighted-Kappa-based loss highly similar models, and thus disagreement between the in- function, which allowed to improve the performance of a DR ferences is indicative of the global model uncertainty (Gal and grading network comparing to training with the cross-entropy Ghahramani, 2016). The authors state that the Bayesian treat- loss. The authors tested di erent input images sizes (squares ment worked better than the softmax output as a proxi for un- with side of 128, 256, 384 and 512 pixels), and the highest res- certainty and also show that uncertainty-aware decision referral olution yielded the best performance. The network, which has can improve the diagnosis. 11.3 million parameters for the 512512 case, was trained in a public dataset (Kaggle DR detection). To overcome the high unbalancing of the dataset, the data is artificially balanced using 1.2. Contributions data augmentation, oversampling the least represented classes. Despite the high classification performance of these systems, Many research work has already addressed DR grading, ex- their black box behavior hinders the application on screening plainability and uncertainty estimation separately. Instead, in settings. Indeed, providing explanations for DR grading is chal- this study we propose DRjGRADUATE, a system that is capa- lenging due to the complexity of the task, and thus research ble of simultaneously providing an explanation and a measure is mainly focused on the binary detection and referral tasks. of uncertainty together with a DR severity grade label. Our For instance, Quellec et al. (2017) analyzed several explana- main contributions are the following: tion methods for the referable DR decisions of their CNNs’. Namely, the authors compare deconvolution (Zeiler and Fergus, estimation of a DR grade with an associated uncertainty, 2014), sensitivity analysis (Simonyan et al., 2014) and layer- thus allowing to better establish which cases require fur- wise relevance propagation Bach et al. (2015). Based on (Si- ther inspection by specialists, which is of high interest in monyan et al., 2014), the authors also proposed the creation of computer-aided diagnosis systems; a pixel-wise explanation map that is jointly optimized with the CNN prediction during training. Costa et al. (2019) proposed  explanation of the algorithm decision, easily and directly the attention map generation via a pre-trained network coupled interpretable by the retinal specialists, in the form of an at- with a Multiple Instance Learning (MIL)-based approach for tention map, which will leverage its use by medical profes- DR detection. Namely, the network is cropped at early layers sionals since it partially mitigates the black-box behaviour to reduce the model’s receptive field, and the classification is associated with deep nets; 4 Teresa Araujo ´ et al. / Medical Image Analysis (2020) extensively validated deep learning-based method for DR 2.1. Predicting a grade with uncertainty severity grading of eye fundus images in the internation- For each input image I, DRjGRADUATE predicts a DR ally recognized DR levels. grade y ˆ and a grade uncertainty u. Let M be the output of the last layer of DRjGRADUATE’s backbone after assessing an input image I. M has size N  N  F, where N and F are, 2. DRjGRADUATE for DR grading respectively, the side and number of features. The lesion map L is computed by performing a 1  1  1 convolution with a The backbone of DRjGRADUATE, depicted in Fig. 2, is linear activation over M, so that L 2 R has size N  N  1. composed of several convolutional-batch normalization blocks L is a map where the value of each element L indicates the i j interleaved with max-pooling layers that has already shown po- y presence of a lesion of grade bL e on a patch of size depen- i j tential for DR grading tasks (de la Torre et al., 2017). The top dent on the model’s receptive field in the input. Following the of the network is composed of a novel component that allows MIL approach, the grade of an image is initially estimated as for the inference of the image’s grade y ˆ together with the cor- y ˆ = max(L) (Costa et al., 2019). responding explainabilityE and the explicit uncertainty u of the Having into account assumptions 1 and 3, y ˆ 7! y ˆ must r g decision, without using labels other than the DR grade, such as ensure that p > p > p > : : : , and p > p > by ˆ c by ˆ c1 by ˆ c2 dy ˆ e dy ˆ e+1 r r r r r lesion-wise annotations. p > : : : , i.e., that the classes near to y ˆ have higher prob- dy ˆ e+2 r Let DRjGRADUATE be a classification model that predicts ability than the others. For each image, DRjGRADUATE must y ˆ = argmax (p), where p is the class-wise probability and g c thus predict a generalized Bernoulli distribution biased around p = 1, with C denoting the number of considered classes. c the classes closer to y ˆ . This bias is introduced in the model by The following assumptions regarding DR grading are made: means of a Gaussian distribution N centered on y ˆ and thus the [assumption 1] grading is a discrete ordinal classification image-wise grade probability p is: problem with an image label y 2 G = f0; 1; 2; 3; 4g, where G is the set of possible DR grades; in this scenario, we want ( ! ) 1 [i y ˆ ] to optimize DRjGRADUATE by minimizing the distance be- r p = exp 0:5 ; 8c 2 G tween the target and the prediction, min(jy y ˆ j), and [assump g 2 tion 2] there is always a di erent lesion type per grade on the p = P (1) universe of possible DR-related lesions, DR ; thus, if one lesions c2G c associates a numerical value to each lesion from DR , ac- lesions cordingly to the DR stage in which they appear, and thus y ˆ where  is the variance of N and p 2 [0 ; 1]. Fig. 3a illus- can be predicted by following the standard MIL assumption, trates the sampling of the probability values p from the Gaus- i.e. max (DR ) 7! y ˆ . lesions g sian probability distribution. Knowing that the uncertainty of DR grading can be seen as a regression problem and con- the prediction can be measured by the value of its entropy (Fris- sequently the predicted image grade y ˆ can be inferred from a ton, 2010), i.e. u is directly related to S = p log(p ), it fol- c c 2 2 continuous output y ˆ 2 R . This also allows for the implicit r + lows that since S depends on  , then u depends on  . Thus Multiple Instance Learning (MIL) of the problem, as we can to by estimating  it is possible to infer the uncertainty associated retrieve y ˆ from a predicted lesion map L, y ˆ = max(L). Further- r r with each predicted grade, y ˆ = argmax(p). more, the uncertainty u could be estimated by u / jby ˆ e y ˆ j . r r In DRjGRADUATE,  is inferred by using a fully- However, this uncertainty does not properly hold in scenarios connected layer on top of M. Since both y ˆ and  de- where the regression cut-o thresholds are adjusted to better pend on M, both grading and the associated uncertainty af- fit the data distribution and malignancy priority, as in many fect the backbone of the model, easing its convergence. In- Kaggle competition solutions . Also, in terms of clinical ap- deed, DRjGRADUATE is trained end-to-end minimizing the plicability, it is of higher interest to provide a class-wise prob- loss function: ability to clinicians since it is easier to interpret. Because of this, [assumption 3] DRjGRADUATE assumes that the pre- L = [b log(p )] + (1 ) (2) c c diction and its uncertainty can be modeled by a Gaussian curve c2G centered near the most probable class. Namely, for each im- where b is the hot-encoded (i.e., a binary representation of) i age, DRjGRADUATE predicts a variance  , which allows for and 2 [0; 1] is a weighting factor. The first term, the categori- y ˆ 2 R 7! y ˆ = argmax(p). With this, DRjGRADUATE is r + g cal cross-entropy, leads DRjGRADUATE to predict y ˆ and thus 2 r capable of providing an uncertainty  and a class-wise prob- y ˆ as close to the truth as possible. The second term regularizes ability that still attempts at minimizing the distance between to avoid excessive Gaussian spreads and thus the stagnation the prediction and the ground truth. The next sections describe of the learning on the local minimum of near equiprobable pre- DRjGRADUATE with more detail. dictions. The global loss L and the structure of DRjGRADUATE im- plicitly introduce an inter-dependency between y ˆ and  that be denotes rounding to the nearest integer, bc rounding down and de allows to deal with images of di erent complexity. For less rounding up. complex images, both  and the preliminary class prediction https://storage.googleapis.com/kaggle-forum-message- attachments/88655/2795/competitionreport.pdf error, jby ˆ e y ˆ j, should be low. However, for dicult images, r r Teresa Araujo ´ et al. / Medical Image Analysis (2020) 5 MP Lesions Conv ... 2×(1×1) N N 1 640 64 FC Uncertainty Maxpool 2×2 Block i Convolution 3×3 + Batch norm + ReLU Convolution 3×3 + Batch norm + ReLU Fig. 2: DRjGRADUATE’s architecture. The network is composed of several convolutional-batch normalization blocks interleaved with max-pooling layers. For each input image I, DRjGRADUATE predicts a diabetic retinopathy (DR) grade y ˆ and a grade uncertainty u. M is the output of the last layer of DRjGRADUATE’s backbone after assessing I. The lesion map L is computed by performing a 1  1  1 convolution with a linear activation over M. L is a map where the value of each element L indicates the presence of a lesion of grade bL e on a patch of size dependent on the model’s receptive field in the input. Following the Multiple i j i j Instance Learning approach, the grade of an image is initially estimated as y ˆ = max(L). The output probabilities (p ; c 2 f0:::4g) are then computed from y ˆ based r c r on a Gaussian distribution centered on y ˆ with variance equal to the uncertainty u, with u learned during training. high values of  can compensate for classification errors since two-dimensional multivariate normal distribution with standard 0 0 the higher spread of the Gaussian tends to increase the predicted deviation  , with  related to the receptive field of the net- probability of the correct class, lowering the cross-entropy loss. work, and s is the stride of DRjGRADUATE. The convolution Fig. 3b shows the evolution of L (Eq. 2) as function of y ˆ and with the Gaussian kernel allows to give higher priority to central , for di erent targets y, y 2 G. These plots show that when pixels of the receptive field, which have a higher contribution to y ˆ is close to y the L loss reaches its minimum. Further, if one the value of each L . r i j freezes y ˆ at a value close to the ground truth, the  that leads 2.3. Dealing with class imbalance to a lowerL loss is smaller than if y ˆ is fixed at a value far from DR grading is an imbalanced problem skewed towards the the ground truth. healthy cases, similarly to many other medical imaging prob- lems. Data imbalance introduces convergence issues during 2.2. Explainability of the DR grade model training, demanding for class balancing strategies. Over- The goal of this study does not include a pixel-level pre- sampling has shown to be an overall good choice for dealing diction of lesions. Even human-performed lesion segmenta- with class imbalance (Buda et al., 2018). However, exces- tion shows high variability, due to the extremely small size of sively reducing class imbalance may lead to an unnecessarily some lesions and also to di erences in the annotation style. In- high trade-o between the model’s sensibility and specificity, deed, some ophthalmologists prefer to only annotate clusters i.e. increase the misclassifications of healthy images, since the of lesions, whereas others are more detailed and prefer to give most represented classes will now tend to be relatively over- lesion-wise markings. Since DR image grading is performed classified. With this in mind, DRjGRADUATE is trained using based on the presence of the typical DR lesions, pinpointing an iterative batch balancing scheme that has already shown suc- relevant lesions should be a valuable tool for doctors to assess cess for DR grading : the algorithm’s prediction without the need for a refined seg- mentation. t1 t1 w = r w + (1 r )w (4) ct co c f DRjGRADUATE provides a grade-wise explanation map, as where w is the expected ratio of images of grade c on the batch ct illustrated in Fig. 4. For each grade c the lesion presence en- at epoch t, w is the representativeness of class c on the dataset, c 0 coded in L is translated to an explainability map E of the same w is the ratio of class c on the batch after f training epochs c f size of I via: and r is the decay rate that controls the update of w . With this, ct in the beginning of the training the network is confronted with > a completely balanced set and does the initial learning based on 1 if c 0:5  L < c + 0:5 i j E (s i; s j) = this, and then progressively starts to adjust (or fine-tune) to a 0 otherwise more real data distribution, but still not as skewed as the origi- E = E ~N 0 (3) c 2 ; nal training set distribution. The oversampling of the least rep- resented classes is performed via online horizontal and vertical where E is a matrix of the same size of E , c 2 G n f0g is one of the possible pathological grades, i; j are each row and column position, respectively, in the L map of side N https://www.kaggle.com/blobs/download/ (i; j 2 N j i; j < N), and N is a convolutional kernel from a forum-message-attachment-files/2797/report.pdf 2 ; Block 1 Block 5 6 Teresa Araujo ´ et al. / Medical Image Analysis (2020) y = 0 y = 1 y = 2 y = 3 y = 4 0.0 0.16 0.31 p 0.47 0.63 0 1 ŷ 2 3 4 Diabetic retinopathy grade 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 (a) Probability Gaussian distribution with ŷ ŷ ŷ ŷ ŷ r r r r r mean y ˆ and standard deviation  (Eq. 1). (b) L loss values (Eq. 2, with = 0:7) as function of y ˆ and , for y 2 G. Loss colorbar: 0 2 . Fig. 3: Uncertainty estimation (Eq. 1) and loss representation (Eq. 2). Messidor-2 Kaggle (test) IDRID SCREEN-DR DMR 0.5<L ≤1.5 Gaussian processing 1.5<L ≤2.5 Gaussian processing N Fig. 5: Class distribution of the DR grading datasets used for evaluation. R0/no DR, R1, R2,  R3, R4/PDR,  SDR,  PPDR. 2.5<L ≤3.5 Gaussian N processing Lesions have been used for evaluation. Fig. 5 shows the class distribu- L>3.5 Gaussian processing tion of these datasets. Kaggle DR detection: This dataset was published by Kag- gle , comprising a large number of high resolution images (ap- Fig. 4: Attention maps construction. Gaussian processing is described in Eq. 3. proximately 35 000 in the training set, 55 000 in the test set). Images were acquired under a large variety of conditions, using di erent types of cameras. Images are labeled by a single clin- flips, rotations, brightness adjustments and contrast normaliza- ician with the respective DR grade, out of 4 severity levels: 1 tion. - mild, 2 - moderate, 3 - severe, and 4 - proliferate DR. Level 0 stands for the absence of DR. The organizers state that this 3. Experiments dataset is similar to any real-world dataset, in which noise is present in both images and labels. Consequently, images may Intra- and inter-datasets experiments were performed for contain artifacts, be out of focus and/or with inadequate expo- evaluating DRjGRADUATE in terms of: 1) DR grading per- sition, aiming at promoting the development of robust and re- formance, 2) uncertainty estimation, and 3) explainability. The liable algorithms that behave well in the presence of noise and DR grading performance was quantitatively evaluated on sev- across a variety of images. eral DR-labeled eye fundus image datasets, and class-wise per- Messidor-2: The Messidor-2 public dataset (Abramo ` et al., formance was analysed. The relation between the uncertainty 2016) consists of color eye fundus images, one fovea-centered and the algorithm’s performance was assessed. Further, the un- image per eye, of 874 subjects with diabetes, for a total of 1748 certainty estimation on good and bad quality images was com- images. Images were acquired using a color video 3CCD cam- pared, using image quality-annotated datasets, in order to assess era on a Topcon TRC NW6 non-mydriatic fundus camera with a whether bad quality images would be associated with higher un- 45 FOV, at 1440960, 22401488, or 23041536 pixels (px). certainties. Finally, uncertainty estimation on unfamiliar data The ground truth (Krause et al., 2018) corresponds to an ad- types was also analysed, aiming at inspecting whether the es- judicated consensus of three specialists. A total of 1744 images timated uncertainties could allow to detect unknown image were used, as the remaining were adjudicated as ungradable. types (i.e., outliers). Di erent medical images were consid- IDRID: The Indian Diabetic Retinopathy Image Dataset ered: colon (colonoscopy), skin, cataracts surgery and breast (IDRiD) dataset (Porwal et al., 2018) is composed of 413 im- biopsy microscopy images. For the explainability, the attention ages of 42882848 px, 50 FOV, acquired using a Kowa VX-10 maps produced by DRjGRADUATE were qualitatively evalu- alpha digital fundus camera. Images were classified accord- ated through visual inspection, as well as quantitatively, by ingly to the International Clinical Diabetic Retinopathy Scale. comparison with lesion annotations performed by an ophthal- mologist. 3.1. Datasets https://www.kaggle.com/c/diabetic-retinopathy-detection https://medicine.uiowa.edu/eye/abramo http://www.adcis.net/en/third- 3.1.1. DR grading datasets party/messidor2/ The Kaggle DR detection training dataset was used for train- https://www.kaggle.com/google-brain/ ing DRjGRADUATE, and several public and private datasets messidor2-dr-grades Probability Teresa Araujo ´ et al. / Medical Image Analysis (2020) 7 Table 2: Modified Davis grading used in Takahashi et al. (2017) for the DMR dataset annotation. Grade Description No DR No signs of diabetic retinopathy Microaneurysm, retinal hemorrhage, hard exudate, SDR retinal edema, and more than 3 small soft exudates Soft exudate, varicose veins, intraretinal microvascular PPDR abnormality, and non-perfusion area over one disc area Neovascularization, pre-retinal hemorrhage, vitreous PDR hemorrhage, fibrovascular proliferative membrane, and tractional retinal detachment Fig. 6: Examples of ophthalmologist’s annotations on the SCREEN-DR dataset.  MA, HEM,  EX, CWS, IrMA,  NV,  PHEM,  PFIB. DMR: The dataset announced in Takahashi et al. (2017) (DMR) is composed of 9939 posterior pole color fundus images DR1: This dataset (Pires et al., 2012) is constituted by 5776 (27202720 px) from 2740 diabetic patients. Images were cap- images with resolution 640480 px, captured using a Topcon tured using a NIDEK AFC-30 fundus camera, and have a 45 TRC-50X mydriatic camera (45 FOV). From these, 1300 are FOV. One to four images may exist from a given patient. Im- labeled as being of good quality (do not contain blur and are ages were graded by modified Davis grading, described in Table centered on the macula) and 1392 have poor quality (blur). 2. The international scale grades (Table 1) were converted to the DRIMDB: The Diabetic Retinopathy Image Database modified Davis scale (Table 2) using the following assumption: (DRIMDB) (Sevik et al., 2014) is composed of 216 images di- NDR = R0, SDR = R1 [ R2, PPDR = R3 and PDR = R4. vided into three classes: good (125), poor (69), and outlier (22). SCREEN-DR: This private dataset consists of retinal images The images were captured by a Canon CF-60UVi fundus cam- from a portuguese DR screening program, managed by the Por- era (FOV of 60 ), and have a resolution of 570760 px. Images tuguese North Health Administration (ARSN). A subset of 966 were annotated by an expert, and good-quality images are suit- images was selected (SCREEN-DR dataset). These were ac- able for use in medical diagnosis by an ophthalmologist. quired using Canon CR-2 AF nonmydriatic retinal cameras and HRF: The High-Resolution Fundus (HRF) database (Kohler ¨ have di erent resolutions (with width and height ranging ap- et al., 2013) contains 18 pairs of images from the same eye proximately from 1500 px to 2600 px). Images from both left from di erent subjects, captured using a Canon CR-1 fundus and right eyes were included, which may be centered at the camera (FOV of 45 ). In each pair there is a poor and a good macula, at the optic disc or elsewhere. Images were graded in quality image. These poor quality images have decreased sharp- 5 levels by a retinal specialist. For 348 of these images, pixel- ness, which can be local or global, due, for instance, to a defo- wise annotations of the lesions MAs, HEMs, CWSs, IRMAs, cused camera. Other quality features, such as image contrast or EXs, NVs, PHEMs and PFIB are available. MAs appear as red- illumination conditions, are not contemplated. dish small and circular dots (Zhang et al., 2010), and may oc- cur alone or in groups (Mahendran and Dhanasekaran, 2015). 3.2. Evaluation HEMs are bright red spots and, contrarily to MAs, their shape 3.2.1. DR grading and appearance has substantial variability (Sil Kar and Maity, 2017). Exudates present a large variety of shapes, sizes and The DR grading performance of DRjGRADUATE is eval- contrast (Harangi and Hajdu, 2013). EXs have a more yellow- uated using the Cohen’s quadratic weighted Kappa (), used ish and bright appearance and sharp edges, whereas CWSs are for measuring inter-rating agreement between raters in ordinal more pale yellow/white without well defined borders (Prentasic multi-class problems. This metric penalizes discrepancies be- and Loncaric, 2016). NVs are thin and weak, exhibiting irregu- tween ratings, in a manner that depends quadratically on the lar appearance such as tortuous patterns and loops (Gupta et al., distance between the prediction and the ground-truth, as fol- lows: 2017). Di erent tools were available for annotating the images, P P C C w O i; j i; j such as dots, circles, polygons and free-hand tool, as exempli- i j = 1 (5) P P C C fied in Fig. 6. The annotated lesions were grouped in the di er- w E i; j i; j i j ent grades to allow a direct correspondence with the predicted where C is the number of classes, w is the weight matrix, attention map for each grade. The R1 ground truth map in- O is the observed matrix and E the expected matrix. w is i; j cludes only the annotated MAs; R2 contains only HEMs; R3 the weight penalization for i; j element, which is given by: contains HEMs, CWSs and IRMAs; R4 includes NVs, PHEMs (i j) th . O is the number of samples from the j class (ground and PFIB. EXs were also annotated but not included in the 2 i; j (C1) th grade maps since they are not directly evaluated in the grading truth) which are classified by the model as being from the i process. class. E is the outer product between the prediction and i; j ground truth classification histogram vectors, normalized such 3.1.2. Quality assessment datasets Three di erent eye fundus images datasets with image qual- ity information were used on uncertainty-related experiments. https://www5.cs.fau.de/research/data/fundus-images/ 8 Teresa Araujo ´ et al. / Medical Image Analysis (2020) P P that E = O .  values range from 1 (complete dis- intensity to a progressively lower value, and evaluating the in- i; j i; j i; j i; j agreement) to 1 (perfect agreement). fluence on the prediction’s uncertainty. However,  alone does not completely reflect the perfor- mance of a DR grading algorithm. When dealing with a very 3.2.3. Explainability unbalanced problem, as DR grading (Fig. 5),  is dominated by With the purpose of evaluating the explanation map, a pre- the most represented classes. Using confusion matrices allows dicted object e corresponds to each connected component from to perform grade-wise evaluation and thus understand better the the grade-wise explanation E after applying a fixed probabil- algorithm’s behaviour. We computed normalized versions of ity threshold. The explanation was quantitatively evaluated as the confusion matrices by dividing each value of the original follows: 1) percentage of all e that overlap with the ground matrix by the sum of the values of each row. truth map of the corresponding grade (O ); 2) percentage of ob j g all e that overlap with the ground truth map of the correspond- MIL assumption validation. Preliminary experiments were ing grade or lower (O ); 3) percentage of correctly predicted ob j performed to test whether the MIL assumption hampered the images for which the maximum activation object overlaps with DR grading performance of the network. This was done by the ground truth map of the corresponding grade (O ); 4) per- max training a network in which the MP layer of DRjGRADUATE centage of correctly predicted images for which at least one e (Fig. 2) is replaced by a fully connected (FC) layer, keeping overlaps with the ground truth map of the corresponding grade the remaining network as it was. No significant performance (O ); 5) percentage of ground truth objects that overlap with class difference was found between the non-MIL and MIL models, the predicted map of the corresponding grade or higher (O ); gt with the first tending to present better R0 detection but a higher 6) percentage of all e that overlap with the ground truth map confusion among the disease cases. Because of this, no further of any grade. (O ). Note that the evaluations 2) and 5) are any experiments for this model were performed. justified by the max-pooling assumption made in our pipeline: it is acceptable to assume that if a given region is detected as being from grade 2, for instance, it may also contain signs from 3.2.2. Uncertainty estimation a lower grade which are not highlighted since only the higher Experiments were performed in order to assess how the pre- grade prevails after the max-pooling. diction performance relates with the uncertainty estimation of those predictions. In order to test whether low uncertainty 3.2.4. Statistical tests images are indeed associated with a higher performance, we defined a set of thresholds on the uncertainty values and, for Statistical tests were performed to assess significant di er- each threshold, select the subset of images with estimated un- ences between the uncertainties inferred for di erent datasets. certainty below that threshold. We then compute the  of our Since these data samples did not follow a normal distribu- model predictions on these subsets. Since very low thresholds tion, which is an assumption on most parametric statistical hy- may lead to the inclusion of a number of images which do not pothesis tests, the non-parametric Kruskal-Wallis H Test was truly represent the dataset, and the  computation for these few used (Kruskal et al., 2012). Further analysing the p-value may samples would not be reliable, the uncertainty thresholds start at not be enough when one is dealing with suciently large sam- 0.15 to ensure an acceptable image subset size. This threshold ples, since a statistical test will indicate a significant di erence was determined considering the dataset with the smaller num- in most of the cases (Sullivan and Feinn, 2012), i.e. very large ber of images (SCREEN-DR), for which an uncertainty thresh- sample sizes may lead to the detection of di erences that are old of 0.15 lead to the inclusion of 5% of the dataset, corre- quite small and possibly non relevant. It is thus important to re- sponding to 48 images. Furthermore, cumulative histograms of port the p-value results but also the e ect size, which measures the class-wise number of images as function of the uncertainty the magnitude of the di erence between groups. The e ect size were produced. This allows to exclude the possibility of varia- is defined as ES =   , where  and  are the averages 1 2 1 2 tions in  being influenced by the different proportion of images of the two groups. The Cohen’s d was used as an e ect size per grade instead of being related with the uncertainty. measure, and can be defined as follows: Finally, matrices of the estimated uncertainty were computed by averaging the DRjGRADUATE’s produced uncertainties for d = ES=S (6) all images from each position of the confusion matrix (predicted where S , the pooled standard deviation, is given by: grade vs ground truth grade). 2 2 (n 1) s + (n 1) s Sensitivity analysis. A sensitivity analysis was performed to 1 2 1 2 S = (7) assess the image quality’s influence in the uncertainty esti- n + n 2 1 2 mation. Namely, the influence of image blur was studied by where n and n are the sample sizes for the two groups and s smoothing the input image with Gaussian filters of increas- 1 2 1 and s are the standard deviations for the two groups. ing standard deviation values and inspecting the effect on DRjGRADUATE’s uncertainty estimation. Likewise, contrast The interpretation of d, considering the thresholds proposed was studied since lower contrasts are also characteristic of bad by Cohen and revised by Sawilowsky (2009), can be established quality images. For that, the image’s contrast was decreased via as follows: 1) d = 0:1: very small 2) d = 0:2: small 3) d = 0:5: linear histogram stretching, by setting the image’s maximum medium 4) d = 0:8: large 5) d = 1:3: very large 6) d = 2: huge. Teresa Araujo ´ et al. / Medical Image Analysis (2020) 9 Having into account that our datasets are commonly on the 4.2. Uncertainty estimation order of more than 1000 samples, the e ect size is reported to- The evolution of  with the uncertainty level threshold for gether with the significance level. Since the significant di er- four di erent datasets is shown in Fig. 8a. For the Kaggle DR ence will most likely be achieved in these datasets, d provides detection test set, a  of approximately 0.82 is achieved for the a measure of the magnitude of these di erences. 50% lowest uncertainty images. Likewise, considering the 75% lowest uncertainty predictions,  is of approximately 0.8. The average uncertainty matrix is shown in Fig. 8c, Fig. 8b depicts 3.3. Parameter setting the grade-wise normalized cumulative histogram as function of the inferred uncertainty. Finally, Fig. 9 shows examples of Kag- The input size to the network (Fig. 2) is 640640 3, and the gle DR detection test set images along with DRjGRADUATE’s side of the output feature map, N, is 20. The receptive field, RF uncertainties of the predictions. of the network is of 156, and the stride, s, is of 32. Regarding 0 The uncertainty estimation results for three quality-labeled the attention maps generation, the  of N 0 from Eq.3 was 2; datasets (section 3.1.2) are shown in Fig. 10a. For all datasets, set to RF/6, since it allows to fit the majority of the Gaussian the inferred uncertainty is statistically higher for the bad quality distribution in the kernel window. The weight in the loss subsets (p-value < 0:05), and the respective Cohen’s d is huge (Eq.2) was set to 0.7, in order to preserve the greater importance for the DRIMDB and medium/large for the other datasets. of the log-loss part but still impose some level of regularization. Fig. 10b shows the uncertainty behavior of DRjGRADUATE Regarding the class balancing, we chose r = 0.99 and w = for datasets of unfamiliar data types along with the eye fun- (0:5; 2; 2; 3; 3) ( f = 300) in Eq. 4, aiming at not overfitting the dus image datasets used for evaluation. Between the unfamiliar original dataset distribution and to increase the generalization data types the statistical hypothesis tests indicate that, for each capability of the network. dataset, the results are statistically di erent (p-value < 0:05) from the neighbouring dataset, and the Cohen’s d values sug- 3.4. Training details gest a large to very large magnitude of the di erences. From the BACH to the SCREEN-DR dataset a statistical di erence is Images were cropped around the FOV and resized to the found, with Cohen’s d indicating a medium magnitude of the input size of the network, 640 640 pixels. The model was di erences. Regarding the eye fundus images, some datasets trained for 240 epochs with the Kaggle DR detection training show a statistically di erent average uncertainty. However, in set, using a batch size of 30, the Adam optimizer, and a learning these cases, Cohen’s d values indicate a small to very small rate of 310 . Experiments were performed on an Intel Core magnitude of the di erences. i7-5960X, 32Gb RAM, 2GTX1080 desktop with Python 3.5, Keras 2.2 and TensorFlow 1.8. 4.2.1. Sensitivity analysis Fig. 11 shows examples of images blurred with the Gaus- sian filters of increasing standard deviation (std). Fig. 12 de- 4. Results picts the results obtained for this test in 5 different datasets (Kaggle DR detection test set, Messidor-2, IDRID, SCREEN- DR and DMR), where the curve represents the average of the 4.1. DR grading uncertainty values for each dataset, and the shaded area corre- sponds to a confidence interval of 95%. The sensitivity analysis DRjGRADUATE achieves similar  to other state-of-the-art via contrast adjustment did not suggest any relation between methods, as depicted in Table 3. Note that the majority of image contrast and uncertainty estimation. the state-of-the-art methods did not test on di erent datasets, and thus an assessment of their generalization capabilities is 4.3. Explainability not possible. DRjGRADUATE was evaluated in di erent test datasets and the consistency of the results suggests that the sys- Several examples of DRjGRADUATE’s attention maps vs tem has properly learned discriminant features for DR grading. the specialist’s annotation on SCREEN-DR dataset images are The grade-wise normalized confusion matrices (in percentage) shown in Fig. 13. In each predicted map, the region with max- are shown in Fig. 7. In general, all the matrices have a diagonal imum activation is surrounded by a square. Figs. 13a, 13c, tendency, i.e., the predictions are close to the ground truth. R4 is 13e and 13f illustrate examples where the maximum value of the most misclassified grade. For instance, for both the IDRID the explanation map properly explains the predicted grade. The and SCREEN-DR datasets (Fig. 7c and 7d), the majority of the remaining images show cases of incorrect predictions and/or R4 images is confounded with R3. From the misclassifications incorrect explanations. reported in Fig. 7, the percentage of over-diagnosis is between Table 4 shows the explanation-wise quantitative performance 24% (Kaggle DR detection test set) and 55% (SCREEN-DR) in terms of the number and type of highlighted lesions, O , ob j g for all datasets, which is higher than the 21% expected for hu- O , O , O , O and O (refer to Section 3.2.3). Fig. 14 ob j max class gt any man observers (Krause et al., 2018). Consequently, the per- shows the values of O and O obtained through the appli- max class centage of under-diagnosis is between 45% (SCREEN-DR) and cation of di erent thresholds to the attention maps. The higher 76% (Kaggle DR detection test set), being always lower than the threshold the smaller the objects on the map, resulting in the 79% reported for human observers. higher specificity and lower sensitivity. 10 Teresa Araujo ´ et al. / Medical Image Analysis (2020) Table 3: Quadratic weighted kappa for DRjGRADUATE’s DR grading in four di erent datasets. / 7025 images from the Kaggle DR detection training dataset, * 7000 images from the Kaggle DR detection training dataset, ? 4 classes,  496 out of the 9443 images. Kaggle test Kaggle Messidor-2 IDRID SCREEN-DR DMR? DRjGRADUATE (proposed) 0.74 0.75/ 0.71 0.84 0.74 0.78 de la Torre et al. (2018) 0.74 - - - - - Gonzalez-Gonzalo ´ et al. (2018) - 0.72* - - - - Takahashi et al. (2017) - - - - - 0.79 Predicted Predicted Predicted Predicted Predicted R0 R1 R2 R3 R4 R0 R1 R2 R3 R4 R0 R1 R2 R3 R4 R0 R1 R2 R3 R4 R0 R1 R2 R3 54 40 81 16 3 0 0 4 1 1 79 16 3 2 0 70 21 7 2 0 78 21 1 0 35 48 17 0 0 16 66 18 0 0 25 40 30 5 0 9 41 41 8 1 22 71 6 1 11 25 49 14 1 2 18 62 18 0 0 6 43 49 2 1 6 45 42 6 5 52 38 5 34 55 4 21 73 6 0 0 10 70 20 0 0 15 72 13 0 14 18 68 2 5 0 0 0 17 40 43 0 0 12 59 29 1 5 30 47 17 2 8 22 33 35 0 (a) Kaggle DR detection test (b) Messidor-2 dataset (c) IDRID dataset (d) SCREEN-DR dataset (e) DMR dataset Fig. 7: Confusion matrixes of DRjGRADUATE’s predictions vs the ground truth DR grades for di erent test datasets. Matrixes were normalized by dividing each value of the original confusion matrix by the sum of the values from the corresponding row (total number of images from that grade in the dataset), and the shown values correspond to percentages. Table 4: Assessment of the produced attention maps (threshold=0.3) in com- Because of this, when inferring images with di erent annota- parison with specialist annotated lesions. These metrics are described in sec- tion criteria or with the presence of artifacts, DRjGRADUATE tion 3.2.3. misclassifies images as can be observed in Fig. 15c and 15d. O O O O O O In other scenarios, misclassifications occur due to failed le- ob j g ob j max class gt any sion detection or similarities between di erent lesion types. 0.506 0.677 0.712 0.784 0.526 0.729 Namely, failures such as predicting a R1 image as R0 can be due to missing a single MA. Indeed, the small size and sub- The threshold 0.3 seems to be a good compromise since it al- tlety of these lesions make their detection particularly chal- lows to retain a large percentage of the lesions while producing lenging (Fig. 16a, where no apparent lesion is seen) and easily confoundable with acquisition or anatomical noise (Fig. 16b). regions not too coarse. Confusions between R1 and R2 are likely related with the di- culty of distinguishing between the red lesions MAs and HEMs 5. Discussion (Fig. 16c and 16d), a task that also challenges medical experts since they can be very similar in size and color. Furthermore, 5.1. DR grading di erent grading standards also account for these classifica- tion errors. Namely, for the DMR dataset, the conversion of DRjGRADUATE’s performance does not degrade between SDR = R1 [ R2 does not holds necessarily true in all cases. the Kaggle DR detection test set, the most similar image- For instance, when an image is R3 by the international scale and grading-wise to the training data, and other datasets, due to the presence of > 20 HEMs, it would fit in the SDR cat- showing the generalization capability of the model. Further- egory of the Davis scale. However, through the conversion, this more, considering the misclassified cases, the 2 higher over- image would fall in the PPDR category, which can justify the diagnosis percentage of DRjGRADUATE in comparison to model’s confusion between R1 and R2 that is verified in Fig. 7e. the expected human performance moderately promotes patient follow-up, challenging specialists to reassess their opinion. De- Classification errors for R2 and R3 grades are mostly due to spite the high performance achieved, DRjGRADUATE tends the grading rules. Indeed, in cases where there are no other to fail to distinguish between R3 and R4 grades. As shown signs of R3 (Table 1), grading depends on HEM counting. in Figs. 15a and 15b, DRjGRADUATE learned to associate However, this task is non-trivial for specialists since HEM may photo-coagulation treatment and laser marks to R4, following appear in clusters, hindering lesion counting and thus introduc- the training set labels. Indeed, these images may not have signs ing noise on the data. Furthermore, the MIL assumption fails of R4, such as neovessels and pre-retinal hemorrhages, and thus for these cases, as there is no novel lesion to find, leading to pre- should either be marked as photocoagulated and removed from diction errors (Fig. 17b). Instead, DRjGRADUATE may have the dataset or be included with the grade corresponding to the associated lesion size to the grade. Alternatively, despite rely- visible lesions. This issue is further exacerbated by the high ing on a single lesion to perform grade inference, the construc- class imbalance of the data (Fig. 5), since the data augmentation tion of the lesion map L depends on the integration of features and class balancing scheme can contribute to the regulation of with di erent receptive fields and thus the lesion with the max- the weights update but not the limited variability of the samples. imum response may indirectly encode lesion counting, serving Ground truth R4 R3 R2 R1 R0 Ground truth R4 R3 R2 R1 R0 Ground truth R4 R3 R2 R1 R0 Ground truth R4 R3 R2 R1 R0 Ground truth R3 R2 R1 R0 Teresa Araujo ´ et al. / Medical Image Analysis (2020) 11 % images and the prediction. This allows the use of the common cat- 4 7 10 27 44 63 83 89 96 97 99 DMR egorical crossentropy as the main term of the loss function SCREEN-DR 5 9 14 28 42 57 73 84 95 97 99 and thus avoids the need for the direct optimization of the , Messidor-2 16 29 43 62 81 89 97 98 99 which requires higher computational power. On the other hand, 6 14 23 41 60 73 86 90 95 96 98 99 Kaggle-test DRjGRADUATE’s performance is inferior to the 0.84  re- 0.925 ported by Krause et al. (2018) on EyePACS-2 dataset (1818 images) (data publicly unavailable at the time of writing). How- 0.900 ever, their model was trained on an expensive dataset composed of over 1.6 million images, and with ground truth labels ad- 0.875 0.875 judicated by multiple ophthalmologists. Further, the authors 0.850 apply an empirically defined cascade of thresholds over their continuous network output in order to define the grade, ensur- 0.825 ing that a baseline sensitivity was achieved for the higher grade classes. Consequently, the system also tends to overgrade im- 0.800 ages, with a normalized over-grading of 72%, which can lead to a high number of unnecessary follow-ups. On the other hand, 0.775 DRjGRADUATE presents a lower over-estimation and thus a 0.750 more moderate trade-o between over and under diagnosis of DR. 0.725 5.2. Uncertainty estimation 0.700 0.15 0.21 0.29 0.4 0.54 0.7 0.89 1.14 Ideally, DRjGRADUATE should infer higher uncertainties uncertainty for misclassified cases. Indeed, high  values occur for low un- (a) Quadratic weighted Kappa vs uncertainty (xx-axis in logarithmic scale) for certainty (Fig. 8a), indicating that the DRjGRADUATE’s uncer- four di erent datasets (DMR: 4 classes). tainty estimation is a valid measure of the quality of the predic- Predicted 1.0 R0 R1 R2 R3 R4 tion. Fig. 8c further indicates that low uncertainty is associated 0.8 with correct or close classification, since the uncertainties are 0.27 0.35 0.39 0.43 0.38 considerably lower in the diagonal, and tend to increase with the 0.6 0.27 0.26 0.25 0.33 0.35 class distance. Furthermore, all grades follow a similar cumula- R0 0.4 0.31 0.31 0.25 0.25 0.28 tive histogram, as shown in Fig. 8b, indicating that the tendency R1 R2 illustrated in Fig. 8a does not result from the di erent propor- 0.38 0.41 0.27 0.20 0.23 0.2 R3 R4 tions of each grade among the analysed image subsets. Low 0.41 0.57 0.38 0.28 0.21 0.0 uncertainty of misclassifications tends to be due to the presence 0.5 1.0 1.5 uncertainty of lesions that greatly resemble another type, characteristic of (b) Cumulative percentage of images from (c) Average predicted uncer- a di erent grade (e.g., large MAs which are similar to HEMs), each grade included per uncertainty threshold, tainty matrix for the Kaggle or to the erroneous identification of acquisition errors as lesions for the Kaggle DR detection test set. DR detection test set. (Fig. 18). The lower  for higher uncertainty cases suggests that Fig. 8: Uncertainty of the predictions in relation with the DRjGRADUATE’s performance. DRjGRADUATE can be used cooperatively in the clinical prac- tice by asking ophthalmologists to review cases that are most likely wrong. Indeed, case selection via uncertainty may sig- as a proxy of the overall relevant lesion density (Fig. 17a). nificantly reduce the workload during the screening pipeline. For instance, discarding the 25% most uncertain images from 5.1.1. Comparison with the state-of-the-art Kaggle DR detection test set allows to achieve a 0.8 , a per- Overall, DRjGRADUATE achieves similar performance to formance similar to human experts (Krause et al., 2018), thus other state-of-the-art methods and shows high generalization allowing DRjGRADUATE to be used as an independent reader capability. Indeed, our model has a  similar to Takahashi capable of identifying dubious cases. et al. (2017)’s method on the DMR dataset (value computed High uncertainty tends to occur in very obscured and blurred by data provided by the authors) while being trained on a dif- images (Fig. 9a, 9b, 9c, 9e, 9f). In these cases, it is not possible ferent dataset and considering a di erent grading scale (inter- to make a diagnosis - at least not fully -, and thus predictions are national instead of Davis). Note that class-wise performance associated with high uncertainty estimations. This will be fur- comparison is limited because confusion matrices are not com- ther discussed in section 5.2.1. Fig. 9d shows an interesting case monly available in the literature. The proposed system achieves of a R4 image correctly predicted by DRjGRADUATE but with a performance similar to de la Torre et al. (2017), while having high uncertainty (0.999). This image shows a very large PFIB only 2% of the parameters and not requiring highly dedicated (PDR sign), occupying almost half of the FOV. However, as re- loss functions. Indeed, DRjGRADUATE’s architecture already ferred, DRjGRADUATE did not correctly learn these R4 signs promotes the distance minimization between the ground-truth (even less of this size) and thus it is not expected to be certain % images/class Quadratic weighted kappa Ground truth R4 R3 R2 R1 R0 12 Teresa Araujo ´ et al. / Medical Image Analysis (2020) (a) y ˆ = 1; y = 2. (b) y ˆ = 0; y = 0. (c) y ˆ = 2; y = 3. (d) y ˆ = 4; y = 4. (e) y ˆ = 2; y = 4. (f) y ˆ = 0; y = 2. g g g g g g (g) y ˆ = 4; y = 0. (h) y ˆ = 2; y = 2. (i) y ˆ = 3; y = 0. (j) y ˆ = 3; y = 4. (k) y ˆ = 1; y = 1. (l) y ˆ = 3; y = 3. g g g g g g (m) y ˆ = 2; y = 2. (n) y ˆ = 0; y = 0. (o) y ˆ = 2; y = 2. (p) y ˆ = 3; y = 2. (q) y ˆ = 3; y = 3. (r) y ˆ = 4; y = 4. g g g g g g Fig. 9: Uncertainty-associated predictions for images from the Kaggle DR detection test set. The color of the frame represents the uncertainty, y ˆ is the prediction and y the ground truth grade. Images were cropped around the field-of-view and resized to a square image. Uncertainty colorbar: 0:05 1:406 in making a diagnosis based on such image. This large white uncertainty (Fig. 9f). structure may have been interpreted as partially obscuring the However, since image quality is not trivial to define within image and thus the high uncertainty. Further, high uncertainties the context of clinical diagnosis, image quality assessment can also result from the resemblance of certain structures in the tends to be a subjective process (Wanderley et al., 2019). One images with DR signs, such as CWSs (Fig. 9i), and with pho- can consider more than two levels of quality, such as good, par- tocoagulation marks (Fig. 9g). Contrastingly, DRjGRADUATE tial, bad, where a partial quality image still has a region suitable correctly predicts, with low uncertainty, the good quality and for diagnosis. However, most datasets are labeled in a binary diagnosable images from Fig. 9k-9r, apart from Fig. 9p. We be- manner (good or bad quality), forcing partially obscured im- lieve this case is an example of the noise present in Kaggle DR ages to considered bad or good. detection dataset annotations, since the image is labeled as R2, Further, as image quality is not the only factor influencing but shows several CWSs (R3 sign), leading DRjGRADUATE to the uncertainty of the model prediction, a full correspondence predict R3 with very low uncertainty (0.05). between image quality and the uncertainty values is not likely. However, a tendency of bad quality images being associated, overall, with higher uncertainties is expected. 5.2.1. Low quality images Poor image quality can lead to the occlusion or distortion DRjGRADUATE produced significantly higher uncertainties of lesions, which hinders a meaningful and complete clinical for bad than for good quality images in each dataset (Fig. 10a), decision (Niemeijer et al., 2006). Lack of focus, sharpness or which suggests it can potentially be used in the DR diagnosis illumination often characterize bad quality images (Pires Dias pipeline without the need for an a priori image quality detec- et al., 2014). In theory, and as suggested by Fig. 9, the algorithm tion step. The HRF datasets is particularly interesting since it should output higher uncertainties for images of poor quality, contains pairs of good/bad quality images, but has the incon- since in the these cases it can not be as sure of the diagnosis venient of having only 18 pairs of images. Some of the HRF as in a good quality image. Indeed, Fig. 9a-f mainly show ex- bad quality images di er only slightly from the good quality amples of largely obscured images which are associated with counterparts (Fig. 19), which leads to similar uncertainty esti- very high uncertainty estimations. Usually, when the image is mation between the two sets. Overall, the bad quality images partially available for diagnosis, if DRjGRADUATE can find seem to still be suitable for diagnosis, which may justify the lesions in the viable image region it predicts a disease grade lower uncertainties for the bad quality subset comparing with with high uncertainty (HEMs from R2, in Fig. 9h), whereas if what is verified on the other datasets. DRIMDB is the more it cannot find any signs it predicts it as healthy, again with high dichotomic dataset, where images are clearly bad or good, ap- Teresa Araujo ´ et al. / Medical Image Analysis (2020) 13 p < 0.05 p < 0.05 d ≈ 0.90 d ≈ 2.79 p < 0.05 d ≈ 0.68 p < 0.05 d ≈ 0.77 p < 0.05 1.2 p < 0.05 d ≈ 0.66 d ≈ 1.48 0.8 p < 0.05 d ≈ 0.85 p < 0.05 0.8 p < 0.05 p > 0.05 p > 0.05 d ≈ 0.09 d ≈ 0.51 d ≈ 0.18 d ≈ 0.05 0.4 0.4 0.0 0.0 Kaggle test DR1 Giana ISIC Cataracts BACH Screen-DR dmr IDRID Messidor-2 HRF DRIMDB Datasets Datasets (a) Uncertainties predicted by DRjGRADUATE for good and bad (b) Uncertainties predicted by DRjGRADUATE for unfamiliar data types from di erent medical quality images from public datasets. p-values obtained through image datasets and for the studied eye fundus images datasets. p-values obtained through the the Kruskal-Wallis H Test and Cohen’s d between each good and Kruskal-Wallis H Test and Cohen’s d between each dataset and its right neighbour on the plot bad sets from a given dataset are presented. are presented. Fig. 10: Uncertainty prediction for di erent types of tasks. mentation scheme performed during DRjGRADUATE’s train- ing, which contemplated brightness adjustments and contrast normalization, and most likely made the network robust to these changes. As the major factor accounting for bad quality in the tests with the bad/good quality datasets (Section 3.1.2 (a) std = 3 (b) std = 6 (c) std = 10 (d) std = 20 and Fig. 10a) was blur/decreased sharpness, not image con- trast/illumination conditions, the sensitivity analysis results are Fig. 11: Examples of convolution with Gaussian filters of increasing standard according to the expectations. deviation (std). 5.2.2. Unfamiliar data types The uncertainty estimation on the unfamiliar data types sug- 1.75 gests that unfamiliar medical image types, i.e., which are not 1.50 eye fundus images, could be detected as outliers by the analysis 1.25 of the uncertainty values produced by the algorithm, which are 1.00 in average higher than for the eye fundus images. 0.75 Colonoscopy images are probably the most similar to the eye 0.50 fundus, showing a similar circular FOV and color (Fig. 10b), 0.25 and sometimes even red and bright structures which may re- 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 semble retinal lesions and anatomical structures. Considering std that the DRjGRADUATE was trained with eye fundus images, Fig. 12: Uncertainty estimation of DRjGRADUATE for images resulting from images that most resemble these can lead to more uncertainty convolution with Gaussian filters of increasing standard deviation (std). The since they are somehow familiar but do not present the usual uncertainty values are the average of the uncertainties for 5 datasets (Kaggle and expected structures for diagnosis, similarly to what happens DR detection test set, Messidor-2, IDRID, SCREEN-DR and DMR) consider- ing a 95% confidence interval. in bad quality eye fundus images. H&E stained microscopy im- ages of breast tissue are the closest to the eye fundus in terms of uncertainty estimation, and are visually far from retinal images. This may be due to the presence of several small rounded struc- parently not presenting a middle term. This is reflected on the tures, darker than the background (cells’ nuclei) which may estimated uncertainties, with the bad quality set showing an av- resemble retinal lesions. Considering that our MIL approach erage uncertainty of 0.750.22 and the good set of 0.300.12. learns to identify the local structures responsible for a given DR Note that both HRF and DR1 contemplate only blur in the bad grade, it makes sense that resemblance at a local level (as hap- quality images, and thus low contrast or bad illumination are pens in microscopy images) leads to lower uncertainties than not considered. resemblance at the global image level (colonoscopy images). Regarding the sensitivity analysis with the blurring opera- tions, we can see based on Fig. 12 that the uncertainty shows a 5.3. Explainability tendency to increase with the increase of the image blur, which decreases image quality. This corroborates our previous find- DRjGRADUATE’s MIL-based architecture leads to the ex- ings that image quality is correlated with uncertainty. When it planation maps to pinpoint regions relevant for diagnosis, and comes to the contrast variation experiments, we believe that the thus highlighting every single lesion is not expected. Because uncertainty did not show any correlation with it due to data aug- of this, the global lesion-wise is not reliable enough to fully Uncertainty uncertainty Uncertainty 14 Teresa Araujo ´ et al. / Medical Image Analysis (2020) Explanation map Annotated lesions Explanation map Annotated lesions MA MA (a) y ˆ = 1; y = 1. (b) y ˆ = 0; y = 1. HEM HEM HEM MA MA (c) y ˆ = 2; y = 2. (d) y ˆ = 3; y = 2. g g CWS HEM IrMA HEM HEM (f) y ˆ = 3; y = 3. (e) y ˆ = 3; y = 3. g NV PHEM NV (g) y ˆ = 4; y = 4. (h) y ˆ = 3; y = 4. Fig. 13: Maps predicted by our method along with the ophthalmologist’s annotated lesions in the SCREEN-DR dataset. The curves are the contours of the explanation maps (threshold=0.3), and the red square indicates the region of most relevance for diagnosis (corresponding to the maximum in the network’s output activation map). Close-ups of relevant regions are shown, with 3 indicating that the region is correctly predicted, and 7 indicating otherwise. Below each close-up the name of the corresponding ground truth lesion (if existent) is shown. Explanation map:  R1,  R2,  R3,  R4; ground truth map:  MA,  HEM,  EX, CWS, IrMA,  NV,  PHEM, PFIB. characterize an image. The detection performance is further that the model detects this lesion type both as a R2 or a R3 hindered by the fact that specialists may not annotate all ab- sign (Fig. 13d), which is understandable since a low number of normalities (Fig. 13) since this can be a large time consuming HEMs is a sign for R2 and a high number a sign of R3. Despite task. On the other hand, DRjGRADUATE is able to correctly this image being annotated as R2, the high number of HEMs identify plausible regions for the reasoning behind the grade of may justify the R3 prediction. However, DRjGRADUATE fails a given image. The analysis of the produced attention maps to detect NVs and other R4 lesions such as PHEMs (Fig. 13g) show the ability of the system to detect lesions as small as MAs and PFIB (Fig. 13h). Most of the times the reason for a R4 pre- (Fig. 13a). As expected, when image quality is poor, image diction is related with the presence of photocoagulation marks diagnosis is hindered and lesions are often missed (Fig. 13b). (Fig. 13g) which corroborates the idea that these were learn by Confusion between MAs and HEMs occurs for some images the algorithm as being R4 signs. Nevertheless, R4 lesions are (Fig. 13c), which is justifiable by the high resemblance of small sometimes highlighted in the maps as signs of disease. For in- HEMs and MAs. The detection of R3 is particularly interest- stance, NVs and PHEMs are detected as R2 signs in these two ing due to the variety of configurations that are associated with images, since they slightly resemble HEMs. this grade (Table 1), which DRjGRADUATE seems to be able Interestingly, the inherent explainability of DRjGRADUATE to identify: a high number of HEMs (Fig. 13d), some CWSs allows for the a posteriori identification of systematic errors (Fig. 13e) or IrMAs (Fig. 13f). Regarding HEMs, one can see during image acquisition. For instance, several Messidor-2 im- Teresa Araujo ´ et al. / Medical Image Analysis (2020) 15 0.80 Omax 0.75 Oclass 0.70 0.65 0.60 0.55 0.50 0.45 0.1 0.2 0.3 0.4 0.5 0.6 0.7 threshold Fig. 14: O and O for di erent threshold values (0.1, 0.3, 0.5 and 0.7) max class (b) y ˆ = 4; y = 4 (Kaggle DR detec- (a) y ˆ = 4; y = 4 (Messidor-2). applied to the attention maps (Eq. 3). tion test). ages have been wrongly classified as R1 due to the presence of dust particles in the camera lens (Fig.18). However, the clear explanation provided would easily allow the clinician to correct the acquisition error and thus improve the overall quality of the screening pipeline. 6. Conclusions (d) y ˆ = 4; y = 0 (Kaggle DR detec- (c) y ˆ = 4; y = 0 (Messidor-2). tion test). Fig. 15: Examples of images from Messidor-2 and Kaggle DR detection test We propose DRjGRADUATE, a novel deep learning-based sets that are predicted by DRjGRADUATE as R4. The curves are the contours approach for DR grading that provides an uncertainty and ex- of the attention maps (threshold=0.3).  R1, R2, R3,  R4. planation associated with each prediction. DRjGRADUATE was trained on the Kaggle DR detection training set, and Acknowledgments achieved state-of-the-art performance in several DR-labeled datasets. The classification of the R4 grade was the least accom- Teresa Araujo ´ is funded by the FCT grant contract plished task, suggesting the training did not correctly capture SFRH/BD/122365/2016. Guilherme Aresta is funded by the the PDR signs. This is most likely due to the type of annotation FCT grant contract SFRH/BD/120435/2016. This work is fi- in the Kaggle DR detection dataset, in which images that under- nanced by the ERDF European Regional Development Fund went photocoagulation treatment and show laser marks appear through the Operational Programme for Competitiveness and to be labeled as PDR even when not presenting the typical signs Internationalisation - COMPETE 2020 Programme, and by Na- of this grade. tional Funds through the FCT - Fundac ¸ao ˜ para a Ciencia ˆ e a The analysis of the predictions’ uncertainties showed the ten- Tecnologia within project CMUP-ERI/TIC/0028/2014. dency of lower  for higher uncertainty cases, thus suggesting that the estimated uncertainty is a viable measure of the quality of the prediction. This means that at test time, i.e. where the References performance of the system is not known, the uncertainty can be Abramo ` , M.D., Lou, Y., Erginay, A., Clarida, W., Amelon, R., Folk, J.C., a measure of how much a specialist can trust the algorithm’s Niemeijer, M., 2016. Improved Automated Detection of Diabetic Retinopa- decision, allowing to rank the predictions and select which im- thy on a Publicly Available Dataset Through Integration of Deep Learn- ages would benefit most from a second evaluation. Further, the ing. Investigative Opthalmology & Visual Science 57, 5200. doi:10.1167/ explanation map produced by our algorithm should be easily iovs.16-19964. interpretable by the retinal specialists, and thus constitutes a Abramo , M.D., Suttorp-Schulten, M.S., 2005. Web-based screening for di- abetic retinopathy in a primary care population: The EyeCheck Project. valuable tool for mitigating the black-box behaviour associated Telemedicine Journal and e-Health 11, 668–674. doi:10.1089/tmj.2005. with deep NNs. 11.668. Acharya, U.R., Ng, E.Y.K., Tan, J.H., Sree, S.V., Ng, K.H., 2012. An inte- Comparing with the state-of-the-art, DRjGRADUATE has grated index for the identification of diabetic retinopathy stages using tex- the advantage of producing both an uncertainty estimation asso- ture parameters. Journal of Medical Systems 36, 2011–2020. doi:10.1007/ ciated with the prediction and an explanation map highlighting s10916-011-9663-8. the most relevant regions for the classification, without compro- Arev ´ alo, J.F., Lasave, A.F., Zeballos, D.G., Bonafonte-Royo, S., 2013. Diabetic Retinopathy, in: Retinal and Choroidal Manifestations of Selected Systemic mising the DR grading performance. We believe that coupling Diseases. Springer New York, New York, NY, pp. 387–416. doi:10.1007/ an uncertainty estimation and explanation with the DR grade 978-1-4614-3646-1{\_}21. will ease the application of DRjGRADUATE in a clinical set- Bach, S., Binder, A., Montavon, G., Klauschen, F., Muller ¨ , K.R., Samek, W., ting. 2015. On pixel-wise explanations for non-linear classifier decisions by layer- hits 16 Teresa Araujo ´ et al. / Medical Image Analysis (2020) (a) y ˆ = 0; y = 1. (b) y ˆ = 1; y = 0. (a) y ˆ = 3; y = 3. (b) y ˆ = 3; y = 2. g g g g Fig. 17: Examples of R2 and R3 Kaggle DR detection dataset images (not used for training) along with DRjGRADUATE predictions. The curves are the contours of the attention maps (threshold=0.3).  R1, R2, R3, R4. (c) y ˆ = 1; y = 2. (d) y ˆ = 1; y = 2. g g Fig. 16: Examples of R0-R2 Kaggle DR detection dataset images (not used for training) along with the DRjGRADUATE predictions. The curves are the contours of the attention maps (threshold=0.3).  R1, R2, R3, R4. (a) y ˆ = 1; y = 0. (b) y ˆ = 1; y = 0. g g Fig. 18: Examples of healthy Messidor-2 images mispredicted as R1. These wise relevance propagation. PLoS ONE 10, 1–46. doi:10.1371/journal. images show camera artifacts (present in the same location across several dif- pone.0130140. ferent images) which are detected as R1 typical lesions. The curves are the Buda, M., Maki, A., Mazurowski, M.A., 2018. A systematic study of the class contours of the attention maps (threshold=0.3).  R1,  R2, R3, R4. imbalance problem in convolutional neural networks. Neural Networks 106, 249–259. doi:10.1016/j.neunet.2018.07.011. Costa, P., Araujo, ´ T., Aresta, G., Galdran, A., Mendonc ¸a, A.M., Smailagic, A., Campilho, A., 2019. EyeWeS : Weakly Supervised Pre-Trained Convolu- Gupta, G., Kulasekaran, S., Ram, K., Joshi, N., Sivaprakasam, M., Gandhi, R., tional Neural Networks for Diabetic Retinopathy Detection, in: International 2017. Local characterization of neovascularization and identification of pro- Conference on Machine Vision Applications (MVA), Tokyo. pp. 5–15. liferative diabetic retinopathy in retinal fundus images. Computerized Med- Cuadros, J., Bresnick, G., 2009. EyePACS: An adaptable telemedicine system ical Imaging and Graphics 55, 124–132. doi:10.1016/j.compmedimag. for diabetic retinopathy screening. Journal of Diabetes Science and Tech- 2016.08.005. nology 3, 509–516. doi:10.1177/193229680900300315. Hall, A.B., 2011. Recognising and managing diabetic retinopathy How I look Dutra Medeiros, M., Mesquita, E., Papoila, A.L., Genro, V., Raposo, J.F., 2015. for diabetic retinopathy : a vision technician s experience. Community Eye First diabetic retinopathy prevalence study in Portugal: RETINODIAB Health 24. StudyEvaluation of the screening programme for Lisbon and Tagus Valley Harangi, B., Hajdu, A., 2013. Improving automatic exudate detection based region. British Journal of Ophthalmology 99, 1328–1333. doi:10.1136/ on the fusion of the results of multiple active contours. Proceedings - Inter- bjophthalmol-2015-306727. national Symposium on Biomedical Imaging , 45–48doi:10.1109/ISBI. Friston, K., 2010. The free-energy principle: A unified brain theory? Nature 2013.6556408. Reviews Neuroscience 11, 127–138. doi:10.1038/nrn2787. Kendall, A., Gal, Y., 2017. What uncertainties do we need in Bayesian deep Gal, Y., Ghahramani, Z., 2016. Dropout as a Bayesian Approximation: Rep- learning for computer vision?, in: Advances in Neural Information Process- resenting Model Uncertainty in Deep Learning, in: Proceedings of the 33rd ing Systems, Long Beach. pp. 5575–5585. International Conference on Machine Learning, New York. doi:10.1109/ Kohler ¨ , T., Budai, A., Kraus, M.F., Odstrcilik, ˇ J., Michelson, G., Hornegger, J., TKDE.2015.2507132. 2013. Automatic no-reference quality assessment for retinal fundus images Garg, S., Davis, R.M., 2009. Diabetic retinopathy screening update. Clinical using vessel segmentation, in: Proceedings of CBMS 2013 - 26th IEEE In- Diabetes 27, 140–145. doi:10.2337/diaclin.27.4.140. ternational Symposium on Computer-Based Medical Systems, pp. 95–100. Gargeya, R., Leng, T., 2017. Automated Identification of Diabetic Retinopa- doi:10.1109/CBMS.2013.6627771. thy Using Deep Learning. Ophthalmology , 1–8doi:10.1016/j.ophtha. Krause, J., Gulshan, V., Rahimy, E., Karth, P., Widner, K., Corrado, G.S., Peng, 2017.02.008. L., Webster, D.R., 2018. Grader Variability and the Importance of Reference Gonzalez-Gonzalo, ´ C., Liefers, B., Ginneken, B.V., Sanchez, ´ C.I., 2018. Im- Standards for Evaluating Machine Learning Models for Diabetic Retinopa- proving weakly-supervised lesion localization with iterative saliency map thy. Ophthalmology 125, 1264–1272. doi:10.1016/j.ophtha.2018.01. refinement , 4–6. 034. Gulshan, V., Peng, L., Coram, M., Stumpe, M.C., Wu, D., Narayanaswamy, Kruskal, W.H., Wallis, W.A., Kruskal, W.H., Wallis, W.A., 2012. Journal of the A., Venugopalan, S., Widner, K., Madams, T., Cuadros, J., Kim, R., American Statistical Analysis. Journal of the American Statistical Associa- Raman, R., Nelson, P.C., Mega, J.L., Webster, D.R., 2016. Develop- tion 47, 583–621. doi:10.2307/2280779. ment and Validation of a Deep Learning Algorithm for Detection of Di- Leibig, C., Allken, V., Ayhan, M.S., Berens, P., Wahl, S., 2017. Leveraging abetic Retinopathy in Retinal Fundus Photographs. Jama 304, 649–656. uncertainty information from deep neural networks for disease detection. doi:10.1001/jama.2016.17216. Scientific Reports 7, 1–14. doi:10.1038/s41598-017-17876-z. Teresa Araujo ´ et al. / Medical Image Analysis (2020) 17 jgme-d-12-00156.1. Takahashi, H., Tampo, H., Arai, Y., Inoue, Y., Kawashima, H., 2017. Applying artificial intelligence to disease staging: Deep learning for improved staging of diabetic retinopathy. Plos One 12, e0179790–e0179790. doi:10.1371/ journal.pone.0179790. The Royal College of Ophthalmologists, 2012. Diabetic Retinopathy Guide- lines. Technical Report December. de la Torre, J., Puig, D., Valls, A., 2018. Weighted kappa loss function for multi- class classification of ordinal data in deep learning. Pattern Recognition Letters 105, 144–154. doi:10.1016/j.patrec.2017.05.018. de la Torre, J., Valls, A., Puig, D., 2017. A Deep Learning Interpretable Classi- fier for Diabetic Retinopathy Disease Grading. arXiv . Wanderley, D.S., Araujo, T., Carvalho, C.B., Maia, C., Penas, S., Carneiro, A., (a) Bad quality image, u=0.415. (b) Good quality image, u=0.295. Mendonca, A.M., Campilho, A., 2019. Analysis of the performance of spe- cialists and an automatic algorithm in retinal image quality assessment, in: 2019 IEEE 6th Portuguese Meeting on Bioengineering (ENBENG), IEEE. Fig. 19: Example of a pair of good/bad quality images from the HRF dataset. doi:10.1109/ENBENG.2019.8692506. Wild, 2004. Estimates for the year 2000 and projections for 2030. World Health 27, 1047–1053. doi:10.2337/diacare.27.5. Mahendran, G., Dhanasekaran, R., 2015. Investigation of the severity level 1047DiabetesCareMay2004vol.27no.51047-1053. of diabetic retinopathy using supervised classifier algorithms. Computers Wu, L., Fernandez-Loaiza, P., Sauma, J., Hernandez-Bogantes, E., Masis, M., and Electrical Engineering 45, 312–323. doi:10.1016/j.compeleceng. 2013. Classification of diabetic retinopathy and diabetic macular edema. 2015.01.013. World journal of diabetes 4, 290–4. doi:10.4239/wjd.v4.i6.290. Massin, P., Chabouis, A., Erginay, A., Viens-Bitker, C., Lecleire-Collet, A., Zeiler, M.D., Fergus, R., 2014. Visualizing and Understanding Convolutional Meas, T., Guillausseau, P.J., Choupot, G., Andre, ´ B., Denormandie, P., Networks, in: Fleet, D., Pajdla, T., Schiele, B., Tuytelaars, T. (Eds.), Com- 2008. OPHDIAT c : A telemedical network screening system for diabetic puter Vision – ECCV 2014, Springer International Publishing, Cham. pp. retinopathy in the Ile-de-France. Diabetes and Metabolism 34, 227–234. 818–833. doi:10.1016/j.diabet.2007.12.006. Zhang, B., Wu, X., You, J., Li, Q., Karray, F., 2010. Detection of microa- Nguyen, T.T., Wong, T.Y., 2009. Retinal vascular changes and diabetic neurysms using multi-scale correlation coecients. Pattern Recognition 43, retinopathy. Current Diabetes Reports 9, 277–283. doi:10.1007/ 2237–2248. s11892-009-0043-4. Niemeijer, M., Abramo ` , M.D., van Ginneken, B., 2006. Image structure clustering for image quality verification of color retina images in diabetic retinopathy screening. Medical Image Analysis 10, 888–898. doi:10.1016/ j.media.2006.09.006. Pires, R., Jelinek, H.F., Wainer, J., Rocha, A., 2012. Retinal image quality analysis for automatic diabetic retinopathy detection. Brazilian Sympo- sium of Computer Graphic and Image Processing , 229–236doi:10.1109/ SIBGRAPI.2012.39. Pires Dias, J.M., Oliveira, C.M., Da Silva Cruz, L.A., 2014. Retinal image qual- ity assessment using generic image quality indicators. Information Fusion 19, 73–90. doi:10.1016/j.inffus.2012.08.001. Porwal, P., Pachade, S., Kamble, R., Kokare, M., Deshmukh, G., Sahasrabud- dhe, V., Meriaudeau, F., 2018. Indian Diabetic Retinopathy Image Dataset (IDRiD). IEEE Dataport . Prentasic, P., Loncaric, S., 2016. Detection of exudates in fundus pho- tographs using deep neural networks and anatomical landmark detection fusion. computer methods and programs in biomedicine 137, 281–292. doi:10.1109/ISPA.2015.7306056. Quellec, G., Charriere, ` K., Boudi, Y., Cochener, B., Lamard, M., 2017. Deep image mining for diabetic retinopathy screening. Medical Image Analysis 39, 178–193. doi:10.1016/j.media.2017.04.012. Sawilowsky, S.S., 2009. New E ect Size Rules of Thumb. Journal of Modern Applied Statistical Methods 8, 597–599. doi:10.22237/jmasm/ Schully, T., 2012. Diabetes in Numbers. Nature 485, S2–S3. doi:10.1038/ nn.2369. Scotland, G.S., McNamee, P., Philip, S., Fleming, A.D., Goatman, K.A., Prescott, G.J., Fonseca, S., Sharp, P.F., Olson, J.A., 2007. Cost-e ectiveness of implementing automated grading within the national screening pro- gramme for diabetic retinopathy in Scotland. The British journal of oph- thalmology 91, 1518–23. doi:10.1136/bjo.2007.120972. Sevik, U., Kose, ¨ C., Berber, T., Erdol, ¨ H., 2014. Identification of suitable fundus images using automated quality assessment methods. Journal of Biomedical Optics 19, 046006. doi:10.1117/1.jbo.19.4.046006. Sil Kar, S., Maity, S., 2017. Automatic Detection of Retinal Lesions for Screen- ing of Diabetic Retinopathy. IEEE Transactions on Biomedical Engineering 9294, 1–1. doi:10.1109/TBME.2017.2707578. Simonyan, K., Vedaldi, A., Zisserman, A., 2014. Deep Inside Convolutional Networks: Visualising Image Classification Models and Saliency Maps, in: International Conference on Learning Representations (ICRL), Calgary. Sullivan, G.M., Feinn, R., 2012. Using E ect Sizeor Why the P Value Is Not Enough . Journal of Graduate Medical Education 4, 279–282. doi:10.4300/

Journal

Electrical Engineering and Systems SciencearXiv (Cornell University)

Published: Oct 25, 2019

There are no references for this article.