Red-lesion extraction in retinal fundus images by directional intensity changes’ analysis

Diabetic retinopathy (DR) is an important retinal disease threatening people with the long diabetic history. Blood leakage in retina leads to the formation of red lesions in retina the analysis of which is helpful in the determination of severity of disease. In this paper, a novel red-lesion extraction method is proposed. The new method firstly determines the boundary pixels of blood vessel and red lesions. Then, it determines the distinguishing features of boundary pixels of red-lesions to discriminate them from other boundary pixels. The main point utilized here is that a red lesion can be observed as significant intensity changes in almost all directions in the fundus image. This can be feasible through considering special neighborhood windows around the extracted boundary pixels. The performance of the proposed method has been evaluated for three different datasets including Diaretdb0, Diaretdb1 and Kaggle datasets. It is shown that the method is capable of providing the values of 0.87 and 0.88 for sensitivity and specificity of Diaretdb1, 0.89 and 0.9 for sensitivity and specificity of Diaretdb0, 0.82 and 0.9 for sensitivity and specificity of Kaggle. Also, the proposed method has a time-efficient performance in the red-lesion extraction process.

Diabetic retinopathy (DR), a common disease among the patients with the history of diabetes, is the most important reason for blindness in developed countries. Studies show that approximately 80 percent of patients with 10 years' history of diabetes are influenced with DR. Currently, 93 million people live with diabetes around the world. The early diagnosis and monitoring of DR can prevent from irreversible damages to human eye. DR is defined as the damage to blood vessels in retina. In fact, the high blood sugar in the diabetic patients can close the tiny blood vessels of retina and lead them to have leakage. In some cases, the eye makes new blood vessels that may not work well 1 .
Fundus imaging is defined as the process by which a 2 dimensional representation of the 3 dimensional tissue is achieved with the help of reflected eye. With respect to the significant increment in the prevalence of DR especially between aging population, the necessity of retinal examination becomes more important  . Fundus images are widely utilized for evaluating and grading of DR. Retinal fundus imaging can produce beneficial information about diseases such as DR, macular edema and glaucoma 25 . The symptoms of DR in fundus images are manifested as red lesions including Hemorrhages (HMs) and Micro-aneurysms (MAs) and bright lesions like exudates. The manual verification of such abnormalities in the retinal fundus images is a tedious, time-consuming and error-prone task. Therefore, the design of automatic Computer-Aided Diagnosis methods for the rapid and accurate monitoring of DR is of considerable importance.
In order to evaluate the severity of DR, it is necessary to extract the related retinal abnormalities. However, this process encounters challenges due to the low quality of retinal fundus images, different shapes and sizes for abnormal lesions, similarity between some parts of blood vessels and red-lesions. Many research works focused on suggesting solutions for the mentioned problems. However, long processing time and unnecessary complicated computations have reduced their efficiency. Thus, to propose new methods with high accuracy, simplicity and speed in the detection of abnormalities of retinal fundus images is considerably beneficial to help ophthalmologist in the evaluation of DR.
In this paper, a novel red-lesion detection method for retinal fundus images is proposed which provides efficiency in terms of speed and accuracy. This method works based on the detection of all boundary pixels with the  38 a method for automated MA detection in retinal fundus images is proposed which consists of pre-processing, candidate extraction, feature extraction and classification phases. In the candidate extraction phase, candidates are extracted using peak detection and region growing. After extracting candidate pixels, region growing is performed to grow candidate pixel into the original MA shape. The features extracted in the next phase include local features like shape and intensity features and profile features. However, a large volume of computations should be performed in this method. In Ref. 39 a method is proposed for the detection of MAs that its main phases are similar to Ref. 38 . Candidate pixels are chosen from pre-processed image and MAs are detected through the intensity analysis along the lines segments centered at the candidate pixels. A method for extracting red-lesions from retinal fundus images is suggested in Ref. 40 . After pre-processing, the optic disc is detected and discarded from the image to not include any red lesion candidate. Then, all the red-lesion candidates are detected based on their contrast and intensity features. Also, dynamic shape features are extracted for the candidates. Then, a classifier is utilized to label the extracted candidates. A method for detecting and classifying abnormalities in retinal fundus images is presented in Ref. 41 . In this method a blob-ness measure is proposed and several intensity and shape features are utilized to identify DR related changes. Then, a Support Vector Machine (SVM) is used to classify the abnormal changes of retina.
With respect to the existing datasets for retinal fundus images, it is possible to mention MESSIDOR, DRIVE, STARE, ROC, DIARETDB0, e-ophtha and DIARETDB1 datasets which are publicly available. Also, the accuracy of the red-lesion detection methods is mostly evaluated with sensitivity and specificity parameters. The sensitivity of a method is determined by its capability to correctly identify lesions and specificity measures the degree of correctly identifying non-lesions. Supplementary Table S1 summarizes the datasets and experimental outcomes of the related existing methods.
The main contributions of our proposed method are summarized in several points.
• In this paper, all the boundary pixels are extracted using simple pixel-wise computations. Then, these boundary pixels are utilized in the red-lesion detection process. • It is possible to control the number of pixels considered as the boundary ones to prevent from producing superfluous pixels. As the number of extracted boundary pixels increases, the required processing time for detecting red-lesion increases. However, if the number of extracted boundary pixels is low, the accuracy of detecting red-lesions is reduced and some red-lesions are missed. Therefore, the number of extracted boundary pixels should be equal to a value which makes a trade-off between processing time and detection accuracy. • Red-lesions are considered as intensity changes in all directions. In fact, an important distinguishing feature of blood vessels is that there is at least one direction in the boundaries of blood vessels along which no significant intensity change is observed. This is while in red-lesions it is possible to observe intensity change along almost all directions. • Therefore, the neighbour windows around the boundary pixels of red-lesions have different intensity-based features. • In this paper, a novel optic disc localization algorithm is proposed to identify the location of optic disc with simple computations. The importance of proposing this algorithm is to remove superfluous boundary pixels which are mistakenly considered as the boundary pixels for red-lesions. With respect to the optic disc zone, it should be noted that no red-lesion is usually observed inside this zone. Therefore, if the optic disc zone is correctly localized, any boundary pixel which is located inside it is removed and the production of false positive cases inside the optic disc zone is prevented.
The experimental results show that the proposed method has an improved performance in terms of speed and accuracy in comparison with the related state-of-the-art research works.

Method
In this section, the proposed method for red-lesion extraction from fundus images is explained in details. The study was approved by the Ethics Committee of National Institute for Medical Research Development (IR. NIMAD.REC.1398.064).
The novel method called Red-lesion Extraction with Directional Intensity Changes' Analysis (REDICA) consists of four main phases. The main phases are pre-processing, boundary pixel determination, red-lesion extraction, and post-processing. Figure 1 presents the block diagram for REDICA method. All phases are explained in the following.
Pre-processing. In this phase, the quality of fundus image is improved and the image is prepared for processing in the next phases. The method used here for contrast enhancement is Contrast Limited Adaptive Histogram Equalization (CLAHE). It should be mentioned that CLAHE originates from Adaptive Histogram Equalization (AHE) which in turn is a kind of histogram equalization. It should be mentioned that AHE utilizes several histograms each of them is related to a distinct region of image to improve the contrast of image. Using AHE method, the local contrast of image is improved and the edges in each region of image are enhanced. However, the disadvantage of AHE method is to over-amplify the noise in the process of contrast enhancement. This issue is solved in CLAHE method by the limited amplification of contrast. This method has been utilized in a large number of research works with similar purposes [45][46][47][48] . At first, the green channel of the fundus image is extracted and then CLAHE method is applied on it. The reason for extracting green channel is that among red, blue and green channels of image, green channel has the best contrast. Since the proposed method works based on the analysis of intensity, the provision of improved contrast for image is very important and significantly helps in www.nature.com/scientificreports/ the process of red-lesion extraction. Between red and green channels, green channel has a significantly better contrast and the improved contrast of image leads to the provision of more boundary pixels in the boundary pixel determination phase of the proposed method. This in turn leads to the extraction of red-lesions with higher accuracy. In addition, with respect to the blue channel, it should be mentioned that the blue channel of fundus image is very noisy and it is necessary to reduce the noise level before processing the image for the purpose of red-lesion detection. This noise reduction should be performed through some de-noising method which in turn inflicts an extra computational load.

Boundary pixel determination.
In this phase, all the boundary pixels in the fundus images are extracted.
These pixels include the boundary pixels for blood vessels and red lesions.
In order to explain the method proposed for this phase, several notations should be defined. Let I denote a fundus image with m rows and n columns. The pixel located at (i,j) is denoted with p i,j . The intensity value of p i,j is denoted by X i,j . The boundary pixels are the ones where a significant transition occurs in the intensity values of image pixels. Therefore, it sounds reasonable to evaluate the intensity changes in several different directions and pick the highest intensity changes.
In order to determine the highest intensity changes, we define a neighborhood window for each pixel. Let N i,j denote the neighborhood window for p i,j . This window is considered as a s w × s w diamond centered at p i,j . The average intensity value of neighbors of p i,j is denoted with ave i,j . It can be computed via the following equation.
In order to verify the directional intensity changes of the image, we consider two directions the first of which is 135 degrees and is called Direction 1. Also, the second one is 45 degrees and is called Direction 2. In order to verify image in Direction 1, for each pixel we compute the difference between the average intensity of two windows consecutively located at the north-west and south-east sides of the pixel. Let Diff 1 i,j denote the absolute difference between the average values of two near windows around p i,j in Direction 1. In fact, the following equation is true.
It should be mentioned that the pixels which maximize the value of Diff 1 i,j are the boundary pixels in Direction 1. In fact, p i,j is considered to be a boundary pixel in Direction 1, if are true. The set of all boundary pixels locally maximizing Diff 1 i,j are denoted with P 1 . In addition, for each pixel we compute the difference between average intensity of two windows consecutively located at the north-east and south-west sides of Direction 2. Let Diff 2 i,j denote the absolute difference between the average values of two near windows around p i,j in Direction 2. In fact, the following equation is true.
It is worth pointing out that the pixels which maximize the value of Diff 2 i,j are the boundary pixels in Direction 2. In fact, p i,j is considered to be a boundary pixel in Direction 2, if are true. The set of all boundary pixels locally maximizing Diff 2 i,j are denoted with P 2 . It is necessary to compute the values of Diff 1 i,j and Diff 2 i,j for all p i,j s to determine whether or not they are boundary pixels for Directions 1 and 2, respectively. The pseudo-code for boundary pixel determination phase is shown in Algorithm 1. In order to restrict the number of boundary pixels and remove extra points corresponding to negligible intensity transitions in the fundus image, we consider a threshold value denoted by th Diff . Then, only the boundary pixels for which the values of Diff 1 i,j or Diff 2 i,j are higher than th Diff , are selected as the boundary ones for Direction 1 or Direction 2, respectively. For simplicity, the conditions of i,j and Diff 2 i,j are computed for p i,j s (line 2 of Algorithm 1). A loop is written in which all pixels are verified to determine whether or not they can be considered as boundary pixels for Direction 1 and Direction 2 (lines 3 to 13 of Algorithm 1). If all the conditions of R 1 , S 1 , T 1 , R 2 , S 2 and T 2 are true, p i,j can be considered as a boundary pixel for both directions (lines 5 and 6). If only the conditions of R 1 , S 1 and T 1 are true, p i,j can be considered as a boundary pixel for Direction 1 (lines 7 and 8). Also, if only the conditions of R 2 , S 2 and T 2 are true, p i,j can be considered as a boundary pixel for Direction 2 (lines 9 and 10).

Red-lesion extraction.
In this section, the unique features of the boundary pixels of red lesions are extracted from the fundus retinal image and are utilized for discriminating them from other boundary pixels.
In the previous section, all the boundary pixels were found through the pseudo-code of Algorithm 1. However, it is necessary to look for the distinguishing properties of red-lesion boundaries to discriminate them from other boundary pixels specially the boundary pixels of blood vessels.
With the accurate analysis of red-lesions in fundus retinal images, it is observed that red-lesions can be seen as intensity changes in almost all directions. However, for blood vessels, there is at least one direction that no intensity change is observed along that direction. Figure 2 presents this point which is considered as the most important distinguishing feature between boundary pixels. In this figure, intensity changes in different directions are shown with different colors. In part a, a red-lesion is presented and it can be observed that in different directions, there are significant intensity changes. However, in part b, a blood vessel is shown and it can be seen that in one direction (yellow direction) no significant change is observed.
This feature can be used in the discrimination between blood vessel boundary pixels and red-lesion boundary pixels. In continue, we explain how this feature is utilized in red-lesion extraction.
Let P denote the set of all boundary pixels obtained in the previous phase using the pseudo-code of Algorithm. 1. All pixels which belong to P should be verified to determine whether they are boundary ones for blood vessels or red lesions. For each p i,j ǫP , several directions are considered and along each direction several neighborhood window are considered to measure the intensity changes in the direction. Let G j 0 :j 1 i 0 :i 1 = g(i 0 : i 1 , j 0 : j 1 ), 1 ≤ i 0 , i 1 ≤ m, 1 ≤ j 0 , j 1 ≤ n denote a part of image with rectangle shape which includes rows from i 0 to i 1 and columns from j 0 to j 1 . At first, consider Direction 0 which makes 0-degree angle with x axis. Let www.nature.com/scientificreports/ W k,0 i,j (k = 0, ±1, ±2, . . . , ±num) denote the kth neighborhood window of p i,j in Direction 0. Let d 0 denote the width of neighborhood window. The neighborhood windows are considered to be square. Also, num denotes the number of neighborhood windows considered at one side of p i,j . W k,0 i,j is determined using the following equation.
where in (4) x 0 denotes the number of first row in W 0,0 i,j . Also, r is a parameter which determines the overlap between adjacent windows. Note that W 0,0 i,j is determined according to the following rule.
In (5) mean is the averaging operator. For instance, mean G j:j+d 0 i:i+d 0 denotes the average of intensity levels of pixels (in g channel) which belong to G j:j+d 0 i:i+d 0 . Also, in (5) min_all denotes the minimum value of mean G and mean G j:j+d 0 i:i+d 0 . After describing neighborhood windows for Direction 0, we consider Direction θ which makes θ degrees angle with x axis. Let W k,θ i,j (k = 0, ±1, ±2, . . . ) denote the kth neighborhood window of p i,j in Direction θ. Note that W 0,θ i,j is equal to W 0,0 i,j in (4). Let x θ and y θ denote the number of first row and the number of first column in W 0,θ i,j , respectively. W k,θ i,j can be obtained through the following equation.
where in (6) and (7) X θ = x θ − kr , Y + θ = y θ + [kr tanθ] and Y − θ = y θ − [kr tanθ] are used for summarizing the equations. www.nature.com/scientificreports/ The pseudo-code for red-lesion extraction phase is presented in Algorithm 2. It should be mentioned that it is not necessary to verify all directions for each boundary pixel. In fact, it is sufficient to verify several directions in the first quarter ( 0 < θ < 90 ) or in the second quarter (90 < θ < 180 ) for each boundary pixel. The reason is that in Section III. B the boundary pixels in Direction 1 (135 degrees) and Direction 2 (45 degrees) have been determined. In other words, for a boundary pixel p i,j maximizing Diff 1 i,j , it is sufficient to verify several directions in the first quarter. Also, for a boundary pixel p i,j maximizing Diff 2 i,j , it is sufficient to verify several directions in the second quarter. Let S Df and S Ds denote the set of all angles considered for verifying different directions in the first and second quarter, respectively. The mean intensity value of W k,θ i,j is denoted by m k,θ i,j . At first, all the required parameters such as r, d 0 , num, S Df and S Ds are initialized. Then, it is necessary to verify all d 0 × d 0 squares in the image from darkness or lightness point of view. For each d 0 × d 0 square, a new parameter is calculated which characterizes the percentage of dark pixels. By "dark pixel", we mean a pixel the intensity of which is lower than a pre-defined threshold. Let us denote such a pre-defined threshold for defining a dark pixel with T dp . A square window is considered as a dark one if the percentage of its dark pixels is higher than a pre-defined threshold denoted by T D . Also, a square window is considered as a light one if the percentage of dark pixels is lower than another pre-defined threshold denoted by T L . Let Pr d 0 a,b denote the percentage of dark pixels in a d 0 × d 0 square window starting from row a and column b. In line 2, we compute Pr d 0 a,b for all d 0 × d 0 square windows and determine whether the related square window is dark or light after comparing Pr d 0 a,b with T L and T D . In line 1, T dp , T L and T D are also initialized.
Then, a loop is executed in which for every pixel p i,j ǫP different neighborhood windows ( W k,θ i,j s) are built (lines 3 to 15 of Algorithm 2) depending on whether the pixel belongs to P 1 or P 2 . It should be reminded that the set of all boundary pixels locally maximizing Diff 1 i,j and Diff 2 i,j are denoted with P 1 and P 2 , respectively. Also, is the set of all boundary pixels obtained in the previous phase using the pseudo-code of Algorithm. 1. If the pixel belongs to P 1 (P 2 ) and does not belong to P 2 (P 1 ), all the windows ( W k,θ i,j s) in the directions of S Df (S Ds ) should be verified (lines 4 to 9 in Algorithm 2). If the pixel belongs to both P 1 and P 2 , all the windows ( W k,θ i,j s) in the directions of S Df or S Ds should be verified (lines 10 to 12 in Algorithm 2). In all cases, it is necessary to verify the windows in Direction 0 ( W k,0 i,j ). In line 13, for each W k,θ i,j , the values of Pr d 0 w i ,w j is found using line 2 and it is determined whether the related W k,θ i,j is dark or light (line 14). If for all directions, W 0,0 i,j is dark and other neighbor windows including W k,θ i,j s (k = 0) are light, p i,j is considered as a boundary pixel for a red-lesion (lines 17 to 19).
Post-processing. In this part, it is necessary to remove pixels which are mistakenly considered as the boundary pixel for red lesions in the previous phase. Here, we suggest a new algorithm for localizing the optic disc in the fundus retinal image. If we can localize the optic disc zone, we can remove the boundary pixels which are located inside this zone. The reason is that there is usually no red-lesion inside the optic disc. The main point utilized in the optic disc localization is that this zone contains the lightest part of image and also the aggregation www.nature.com/scientificreports/ of blood vessels which are comparatively dark. Therefore, if the standard deviation value for the intensity values of the pixels inside optic disc is computed, it has a high value for this zone. It is possible to grid image and perform computations separately for each grid to determine the grid with the highest variance value. Let d od denote the approximate diameter of optic disc. This parameter can be computed from several fundus images. Our purpose is to find a d od × d od square approximately located on the optic disc zone. The pseudo-code for optic disc localization is presented in Algorithm 3.
At first, the required parameters are set. Note that r d denotes the distance between two adjacent d od × d od squares. In fact, r d is equal to the number of rows (columns) between the starting rows (columns) of two adjacent squares. If r d is equal to 1, all the possible squares with d od × d od dimensions in image are considered and tested. However, it may result in huge volume of computations. Therefore, it sounds reasonable to choose a higher value for r d . It is obvious that the number of overlapped rows and columns between two adjacent squares is equal to Then, a loop is executed in which all test-squares are built (lines 2 to 9 in Algorithm 3). By "test-square", we mean a d od × d od square tested for being the optic disc zone. Also, for each test-square the value of standard deviation of intensity values of its pixels is computed (line 5 in Algorithm 3). Let std t denotes the standard deviation of intensity values of a test-square. std t is equal to

Numerical results
In this section, the performance of REDICA method is evaluated. In order to evaluate the performance of the proposed method, a program has been run in MATLAB. Also, the code has been tested for three different datasets including Diaretdb1 49 , Diaretdb0 50 and Kaggle datasets 51 . Table 1 presents the initial values for the required parameters.
In addition, it should be mentioned that the distribution used in CLAHE method is uniform and the clip limit is set to 0.02. Figure 3 presents a fundus image from Diaretdb1 dataset versus its boundary pixels which are obtained through boundary pixel determination phase in Algorithm. 1. In parts (b) and (c), the boundary pixels are obtained for th Diff = 0.12 and th Diff = 0.08, respectively. It is obvious that the more the value of th Diff , the less the number of obtained boundary pixels. Thus, th Diff should be set to a value that provide sufficient number of boundary pixels. The reason is that if the number of boundary pixels is considerably high, the processing time for red-lesion extraction becomes high while it does not improve the accuracy significantly. Figure 4 presents a fundus image and its boundary pixels of red-lesions obtained from the red-lesion extraction phase of REDICA method. As can be seen, the boundary pixels in parts b, c, d are obtained for d 0 = 10, 16, 22, respectively. Also, part e presents all the boundary pixels for red lesions. As obvious from the figure, the boundary pixels for red lesions are well extracted from visual point of view. Figure 5 presents a fundus retinal image versus the locations of red-lesions localized by the proposed algorithm. The localized red-lesions extracted by REDICA algorithm are shown with blue squares. Also, the redlesions localized by the ophthalmologists are shown with purple squares. It can be seen that red-lesions in the REDICA algorithm are extracted with good accuracy. Figure 6 presents a fundus retinal image for which the optic disc zone is localized using the pseudo-code of Algorithm 3. The optic disc zone is shown with a black square in the image.     (8) and (9), TP, TN, FP and FN denote truly classified lesions, correctly classified non-lesions, nonlesions falsely classified as lesions, and lesions falsely classified as non-lesions, respectively. In order to analyze the statistical difference between the results of our proposed method and the existing methods, a threshold is considered for p-value which is equal to 0.05. It should be mentioned that in the datasets used here, the red-lesions are localized and not exactly segmented. Therefore, in order to determine the values of TP, TN, FP and FN, after extracting the boundary pixels related to red-lesions, it is verified whether or not they are included inside a red-lesion localized and annotated by the ophthalmologist in the dataset. If the mentioned boundary pixel is located inside the spot, such spot is considered as a correctly detected red-lesion by the proposed algorithm (TP). If for a red-lesion annotated by the ophthalmologist, no pixel has been extracted by the proposed algorithm, the mentioned spot is considered as a false negative (FN). In addition, if a boundary pixel is not located inside any annotated red-lesion, it is considered as a false positive spot (FP). Furthermore, if a pixel which is not determined as a boundary one for red-lesion is located inside a non-annotated regions of fundus image, it is considered as a true negative (TN). It should be also mentioned that the accuracy parameters including specificity and sensitivity are computed based on the number of lesions.
In continue, we notice to the results related to Diaretdb1 dataset. Diaretdb1 consists of 89 images that 5 images are normal and the others have the sign of DR. The research works chosen for comparison for Diaretdb1 dataset are the research works of Refs. 35,53 . Dark and bright abnormal lesions are detected with the method of Ref. 53 . In the method of Ref. 53 , the whole fundus image is divided into patches and for each patch a series of features are extracted based on texture and morphology of the patch. In order to extract texture feature, rotation-invariant Local Binary Pattern (LBP) is utilized which is a type of LBP texture indicators. The features related to the morphology are extracted based on granulometry technique. Then, several different classifiers including Random Forests (RF), Gaussian processes and SVM utilizes the extracted features to categorize the abnormal pathologies.  www.nature.com/scientificreports/ The capability of method was evaluated on DIARETDB1 in term of detecting hemorrhages and micro-aneurysms.
In Table 2 the results of comparison between REDICA and the methods of Refs. 35,53 are presented. It was shown that the sensitivity and specificity values in the method of Ref. 53 are equal to 0.75 and 0.75, respectively which are significantly lower than those of our proposed algorithm. In addition, the method of Ref. 35 evaluated its performance on DIARETDB1. The sensitivity and specificity values of method of Ref. 35 are equal to 0.88 and 0.91, respectively which are very close to the values of same parameters in our proposed algorithm. It should be mentioned that the required computations in our proposed method are considerably lower compared to the method of Ref. 35 . The reason is that the method of Ref. 35 requires the segmentation of vessels and the decomposition of image into several layers before determining lesions which in turn affects execution speed and time.
The research work chosen for comparison for Kaggle dataset is Ref. 43 . In Ref. 43 a CNN is trained with a number of informative samples from a large number of medical images. A dynamic weight is assigned to each pixel to present its informativeness level. After each round of CNN training, the weight of each pixel is updated. The final CNN is capable to categorize each pixel in a test image. It should be mentioned that the number of images used for training and testing stages in SeS method is equal to 5287 and 1392, respectively. The number of images randomly selected from the same dataset in REDICA is considered to be 1392, too. The values of sensitivity and specificity for REDICA method and SeS method of Ref. 43 are presented in Table 2. As can be observed, specificity and sensitivity values in REDICA method are 0.9 and 0.82, respectively. Also, specificity and sensitivity values in SeS method are equal to 0.85 and 0.83, respectively. It can be seen in Table 2 that the sensitivity values of both methods are very close and the specificity value of REDICA is higher than the SeS method of Ref. 43 .
In the following, we pay attention to the results related to Diaretdb0 dataset. Diaretdb0 includes 130 images 20 of them are normal and 110 of them have the signs of DR. The research works which are selected for comparison for Diaretdb0 dataset are Refs. 33,54 . In Ref. 54 a method is proposed for extracting hemorrhages from retinal fundus images. At first, all dark regions including blood vessels and hemorrhages are extracted from the image. Then, retinal vessels are segmented and gradually removed from the extracted regions. Table 2 presents the values of sensitivity and specificity for REDICA, the method of Ref. 54 and the method of Ref. 33 for Diaretdb0 dataset. It can be observed that the values of sensitivity and specificity for the method of Ref. 54 are equal to 0.95 and 0.82, respectively. Also, REDICA is capable of providing the values of 0.89 and 0.9 for sensitivity and specificity, respectively. The value of specificity in our method is higher than the method of Ref. 54 . With respect to sensitivity value, it should be noted that the method of Ref. 54 requires the exact segmentation of blood vessels which in turn inflicts a large volume of computations while there is no need to the segmentation in REDICA. In Ref. 33 a method for detecting different kinds of abnormalities in retina including red-lesions is proposed. The method called HEM detector, works based on the analysis of several wavelet bands and Hessian multi-scale analysis. Large structures like red-lesions and blood vessels can be detected in high wavelet levels. Also, with the help of Hessian multi-scale analysis, it is possible to find blob-like structures in the image. It can be observed that the values of sensitivity and specificity for REDICA method are better than those of the method of Ref. 33 .
Also, the Receiver Operating Characteristic (ROC) curves of the proposed algorithm are presented for all datasets. These curves can be observed in Supplementary Fig. S1. In order to draw the ROC curves, we have changed the value of T D and verified its effect on the specificity and sensitivity values. The increment of T D makes more difficult conditions for considering a boundary pixel as a boundary pixel for a red-lesion. The reason is that the related patch located in the center should be darker in this case. Therefore, the number of boundary pixels extracted as the boundary pixels for red-lesions is reduced and the SE value also decreases. In such conditions, the number of boundary pixels falsely detected as the boundary pixels for red-lesions also decreases. Thus, the specificity value increases. The ROC curves present such a relationship between the values of SE and (1 − SP). In addition, the Area Under ROC Curve (AUC) is computed for REDICA method for different datasets. It can be observed that AUC values for Diaretdb0, Diaretdb1 and Kaggle datasets are equal to 0.89, 0.87 and 0.9, respectively. In the methods of Refs. 43,53 , AUC values for Diaretdb1 and Kaggle datasets are equal to 0.83 and 0.89, respectively.

Conclusion
In this paper, a novel method for red-lesion extraction from fundus retinal images is proposed. The method works based on the analysis of intensity changes in the image. Firstly, the method improves the contrast and quality of image for processing. Then, the method extracts all the boundary pixels for red lesions and blood vessels through the intensity changes analysis of image in 45 degrees and 135 degrees directions. In order to distinguish boundary pixels of red-lesions from blood vessels, a distinguishing feature for them is utilized. To extract a distinguishing Table 2. SE and SP values of the methods of Ref. 35,53 for Diaretdb1, the methods of Ref. 33,54 , for Diaretdb0, the method 43 for Kaggle, and the REDICA method for all datasets. www.nature.com/scientificreports/ feature for the boundary pixels of red-lesions, the patches around the boundary pixels are considered and the darkness and lightness of such patches is determined. For the boundary pixels of red-lesions, the nearest patch is dark and the other patches are light in all directions. In addition, in the post-processing step, we can localize the optic disc zone and remove the boundary pixels inside the zone. The reason is that there is usually no redlesion inside optic disc zone. The proposed method has been evaluated for three datasets including Diaretdb0, Diaretdb1 and Kaggle. The new method is capable of providing 0.87, 0.89 and 0.82 for the sensitivity in Diaretdb1, Diaretdb0, and Kaggle datasets, respectively. Also, the values of 0.88, 0.9 and 0.9 are provided by the proposed method for the specificity in Diaretdb1, Diaretdb0, and Kaggle datasets, respectively. The evaluations of the proposed method show that it has improved performance compared to the related existing works.