{'fraction_non_alphanumeric': 0.031978841067564316, 'fraction_numerical': 0.012503005530175523, 'mean_word_length': 4.634930410149033, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Retinal Microaneurysms Detection using Local Convergence Index Features

Retinal microaneurysms are the earliest clinical sign of diabetic retinopathy disease. Detection of microaneurysms is crucial for the early diagnosis of diabetic retinopathy and prevention of blindness. In this paper, a novel and reliable method for automatic detection of microaneurysms in retinal images is proposed. In the first stage of the proposed method, several preliminary microaneurysm candidates are extracted using a gradient weighting technique and an iterative thresholding approach. In the next stage, in addition to intensity and shape descriptors, a new set of features based on local convergence index filters is extracted for each candidate. Finally, the collective set of features is fed to a hybrid sampling/boosting classifier to discriminate the MAs from non-MAs candidates. The method is evaluated on images with different resolutions and modalities (RGB and SLO) using five publicly available datasets including the Retinopathy Online Challenge's dataset. The proposed method achieves an average sensitivity score of 0.471 on the ROC dataset outperforming state-of-the-art approaches in an extensive comparison. The experimental results on the other four datasets demonstrate the effectiveness and robustness of the proposed microaneurysms detection method regardless of different image resolutions and modalities.


I. INTRODUCTION
Diabetic retinopathy (DR) is the most common cause of vision loss among people with diabetes and the leading cause of vision impairment and blindness among middle-aged population in the world [1]. Chronically high blood sugar levels from diabetes are associated with increasing damage to the tiny blood vessels in the retina, leading to diabetic retinopathy [2]. DR can cause vessels in the retina to leak fluid or to bleed, and in advanced stages, new abnormal blood vessels may proliferate on the surface of the retina, which can lead to scarring and cell loss in the retina [3]. Small swellings in the retina's tiny blood vessels, called microaneurysms (MAs), occur at this earliest stage of the disease [3], [4]. In digital color fundus images, MAs appear as tiny, reddish isolated dots near tiny blood vessels [4].
The detection and analysis of MAs is considered as one of the most important clinical strategies for the early diagnosis of DR and blindness prevention in a cost-effective health care practice. However, due to the limited number of ophthalmologists and the large number of people that require screening, an automated computer-aided diagnosis tool can significantly improve the efficiency and reduce the costs in a large-scale screening setting [5].
Several methodologies for the detection of MAs have been proposed, most of which were performed in two stages [6]- [25] (see also Table V). In the first stage, several candidates with similar characteristics to MAs are extracted. In the second stage, a set of features is obtained for each candidate and a classification technique is applied for the discrimination of MAs from non-MA candidates. Niemeijer et al. [7] proposed a hybrid approach for the candidate extraction which is a combination of mathematical morphology techniques and a pixel classification method. The true MAs are then detected using a k-Nearest Neighbors (kNN) classifier and a set of shape and intensity features. For the extraction of initial candidates, Mizutani et al. [10] utilized a modified doublering filter in which the average pixel intensity value of the inner ring is compared with that of the outer ring, then the candidates are re-examined by removal of blood vessels. The candidates were classified using a three-layered feed-forward neural network.
In the method introduced by Sanchez et al. [11], a mixture model-based clustering technique is used for the candidate extraction, which is followed by a logistic regression classifier to generate a likelihood for each candidate based on its color, shape, and texture. Zhang et al. [12] used a multi-scale correlation coefficients-based method and a dynamic thresholding technique for candidates extraction. Then a rule-based classification algorithm was employed for accurate detection of MAs. This approach is further improved by including dictionary learning with a sparse representation classifier [15]. In the method by Ram et al. [13] the candidates are selected using a morphology-based approach in which linear structures are extracted in different orientations. In the next step, a successive rejection-based strategy was designed which passes only true MAs while rejecting false classes of clutter.
Giancardo et al. [14] detected MA candidates using a thresholding technique followed by a Radon transformation at various scanning angles. MAs are then located using a support vector machine (SVM) classifier and features extracted from Radon-space. The ensemble-based framework proposed by Antal and Hajdu [17] detects MAs by selecting the optimal combination of different preprocessing techniques and candidate extractors [6], [8], [12]. Lazar and Hajdu [18] proposed a method based on the analysis of locally rotating cross-section profiles, where the local maximum pixels were selected as candidates. This was followed by a peak detection technique applied on 30 produced profiles, and a set of descriptors including the size, height, and shape of the peaks is obtained. The statistical measures of the resulting directional peak descriptors are then used as the set of features for a Bayesian classifier. Zhang [19] introduced new contextual descriptors which are used in combination with intensity and geometric features in a Random Forest classifier to detect MAs. Adal et al. [20] presented a MA candidate selection approach based on scale-invariant interest-points and blob detection theory. They used a local-scale estimation technique and several scale-adapted region descriptors to characterize detected blob regions. The final classification step was performed using a semi-supervised learning approach.
In the method proposed by Shah et al. [21], initial candidates were extracted by removing vessels and local thresholding. A rule-based classifier and a set of statistical features were employed to classify the candidates into MAs and non-MAs. Dai et al. [22] extracted MAs candidates by a vessel removal technique and gradient vector analysis. The candidates are then classified using a class-imbalance classifier and several features including geometry, contrast, intensity, edge, texture and region descriptors. Recently, Seoud et al. [24] introduced a new set of shape features called Dynamic Shape Features which are employed for the detection of red lesions in retinal images. Wu et al. [23] presented a candidate selection approach similar to Lazars method [18] using profile analysis and region growing. For the final classification stage, the authors used several local and profile-based features and different classifiers such as KNN and Adaboost. In the method proposed by Wang et al. [25] candidates were located using a dark object filtering process. Afterwards, singular spectrum analysis is employed to decompose cross-section profiles of extracted objects and reconstruct a new one. A kNN classifier and a set of statistical features of profiles were used to discriminate the MAs from non-MAs candidates.
In addition to the two-stage approaches, Pereira et al. [26] proposed a multi-agent system method for MAs segmentation using gradient patterns and Gaussian fitting parameters in different directions. Quellec et al. [27] modeled the MAs with 2-D rotation-symmetric generalized Gaussian functions and used a supervised template matching technique in waveletsubbands for the MAs detection. The optimal adapted wavelet transform for MA detection was found by applying the lifting scheme framework and then MAs were detected using template matching in the wavelet domain. Although in recent years deep neural networks have gained popularity for the DR detection and in general in the field of computer vision, only a few papers report specifically on the detection of MAs using deep learning techniques [28], [29]. The results of MA candidate extraction and MA detection algorithms presented in the literature are listed and compared in Tables IV and V. Despite the many published approaches, described above [6]- [29], accurate detection of MAs is still a challenging task. The detection of MAs depends on the image properties and the characteristics of the imaging device such as resolution, image modality, inter/intra-image illumination and contrast variation and compression technique. On the other hand, the variations in MA size, shape and proximity to vessels can also increase the difficulty of identifying them in retinal images. In this paper, we propose a novel method for the detection of MAs using local convergence index filters (LCF) and a random undersampling boosting classifier (RUSBoost). The method extracts several candidates and generates a set of features for each candidate depending upon their intensity, shape and LCF responses. The true MA candidates are then selected using a hybrid sampling/boosting classifier to avoid the drawback of imbalanced data learning and to improve the performance of MA detection. The major contributions of this paper can be summarized as follows.
(i) A novel method is proposed for accurate and reliable detection of microaneurysms with the possibility of applying this method in large screening setups. The method outperforms state-of-the-art techniques. (ii) A new MA candidate extraction technique is employed which extracts the potential MA candidates with a notable improvement in sensitivity; (iii) A new set of features based on local convergence index filters is introduced which are then incorporated into a hybrid sampling/boosting classifier; (iv) The proposed approach is evaluated on five datasets with two image modalities. The high performance of our method provides an indication of the detection power of the proposed approach in handling difficult cases such as subtle MAs and the ones close to vessels; (v) We published two new publicly available retinal image databases, RC-RGB-MA (RGB fundus camera) and RC-SLO-MA (scanning laser ophthalmoscope) to evaluate microaneurysms detection methods on different image modalities. Moreover, we provided a software tool for microaneurysm annotation (RC-MAT) helping the experts to collect more labeled MAs. The tool is publicly available as a software package; This paper is organized as follows. In Section II, the proposed method for microaneurysms detection is presented.
In Section III, we show the experimental results with performance evaluations. Finally, we discuss and conclude our method in Section IV.
II. METHODOLOGY Fig. 1 shows the pipeline of the proposed algorithm. After the preprocessing step, the main phases are 1) candidates extraction using multi-scale multi-orientation gradient weighting and iterative thresholding; 2) features extraction using local convergence index filters; and 3) microaneurysm detection by training a RUSBoost classifier.

A. Preprocessing
Since the green channel of retinal images in RGB datasets provides a better contrast between microaneurysm and background, we only use the green channel (I G ) in our experiments. As a result of the acquisition process, very often the retinal images are non-uniformly illuminated and exhibit local luminosity and contrast variability. In order to make the MB detection more robust, each image is preprocessed using the method proposed by Foracchia et al. [30], which normalizes both luminosity and contrast based on a model of the observed image. The local normalization reduces luminosity and contrast variation of retinal images, and it improves the visibility of the lesions. The result of preprocessing on a sample retinal image is demonstrated in Fig. 2.

B. Gradient Weighting
We compute multi-scale multi-orientation weights for each pixel in a normalized image I N based on the gradient magnitude at that pixel. The gradient magnitude values for different scales (σ G ) and different orientation (θ G ) are obtained by the convolution of the normalized image with the rotated first order derivative of the Gaussian kernel: where θ G ∈ {0, π/12, ...π/2} denotes the rotation angle of the Gaussian derivative kernel and the Gaussian kernel is given by G(x, y, σ G ) = (1/ 2πσ 2 G )e (−x 2 +y 2 )/2σ 2 G where σ G ∈ {1, 2, 3, 4, 5}px is the standard deviation of the Gaussian kernel which is used as scaling parameter.
The weight of a pixel is inversely related to the gradient values at the pixel location, so that for pixels with small gradient magnitude (smooth regions), the weight is large, and for pixels with large gradient magnitude (such as on the edges), Orientations =0 /2 = /3 … = /6 … = 0  the output weight is small. The gradient-weighted image is obtained by Samples of gradient-weighted images of a small patch for different scales and orientations are shown in Fig. 3. The final multi-orientation gradient-weighted image I W O is defined as the sum of the individual gradient-weighted images in all orientations for each scale separately via The multi-orientation gradient-weighted images are used in the successive thresholding process to extract candidates. Examples of the obtained results are illustrated in the last row of Fig. 3. To cover all microaneurysms sizes, the multiorientation gradient-weighted images are obtained at several scales. The final multi-orientation and multi-scale gradientweighted image I W OS is obtained by the summation of I W O over all selected scales

C. MA Candidate Extraction
The candidate extraction step is one of the main phases in automated detection of microaneurysms. In this subsection, a suitable candidate selection algorithm is proposed which decreases the false positive rate as well as the complexity and the computation time by reducing the number of objects for further analysis.
In this step, an iterative thresholding process is applied on the multi-orientation gradient-weighted images (I W O ) to obtain a set of binary images I t σ G corresponding to different threshold values (t) in each scale separately. Using the connected component analysis in each binary image, the objects (C t k ) with a hole inside satisfying the area, eccentricity and extent constraints are identified and selected as MA candidates. In our experiments, if the object's area, eccentricity and extent are smaller than 300 px, 0.9 and 0.3, respectively, then the object is included in the set of selected candidates M t σ G . The parameters are set to the values which produce the highest possible sensitivity.
The final set of MA candidates, M, is the union of all selected candidates at different scales and threshold. The detailed process of the MA candidates extraction is given in Algorithm 1. It is worth noting that the proposed candidate extraction technique excludes the points on the vessels, which reduces the detection error in the final classification step. Fig. 5 shows examples of the final extracted candidates.

Algorithm 1 Candidates extraction
where k ∈ {1, ...., Nt} and Nt is total number of components

D. Feature Extraction
Extracting suitable features and descriptors for the candidate regions is an important step for the final classification stage. Since the MAs appear in different colors and sizes, several shape and intensity features are extracted. The feature set is completed by including the responses and the estimated radii of different local convergence index filters (LCF). The rest of this subsection describes the 29 proposed features which we have defined to characterize and classify MAs. Table I contains the description of the features which we use in the proposed system.
1) Intensity-based features: These features are descriptors indicating the darkness of MAs compared to their neighborhood background. The intensity features are extracted at the center of the candidate, inside the whole candidate object and also in a square neighborhood region which is 3 times as large as the candidate area. The neighborhood region is centered on the center of candidate. As shown in Table I for every candidate in the M set, the average, maximum and minimum of the green intensity values are obtained for the candidate and neighborhood regions, separately.
2) Shape-based features: Since MAs are small and they appear as round structures with a diameter less than 125µm, the following shape-based features are extracted for each candidate region: • Area (S Area ): area of candidate region specified by the actual number of pixels (white pixels in Fig. 6a).
• Convex area (S ConA ): area of candidate convex region specified by the actual number of pixels (white and gray pixels in Fig. 6a).
• Solidity (S Sol ): ratio of the area of candidate (S Area ) over the convex area (S ConA ).
• Extent (S Ext ): ratio of S Area to the pixels in the bounding box as shown by red color in Fig. 6b.
• Perimeter (S P er ): distance around the boundary of the region by calculating the distance between each adjoining pair of pixels.
• Circularity (S CirD ): diameter of a circle with the same area as the region which is equal to 4S Area /π.
• Ellipticity (S AxiA , S AxiB ): lengths of the major and minor axes of the ellipse that has the same normalized second central moments as the candidate region. The major and minor axes are depicted by red lines in Fig. 6c.
• Eccentricity (S Ecc ): ratio of distance between the foci (blue stars in Fig. 6c) and the major axis length (S AxiA ) of the ellipse with a same 2 nd moment as the region.
• Euler number (S Eul ): number of objects in the region minus the number of holes in those objects.

3) LCF-based features:
The LCF filters are based on gradient convergence and not intensity and as such can detect low contrast MAs which otherwise would be easily lost in the background noise. Additionally, the convergence evaluation in a regional band allows the reduction of uncertainty caused by noise. We first give explanations about the local convergence index filters and then introduce the set of LCF-based features extracted for each candidate. The convergence index (CI) Notation Name Description Green intensity value at the center of candidate region Gmean, Gmax, G min Mean, maximum and minimum green intensity values of candidate region f 5−7 G N mean , G N max , G N min Mean, maximum and minimum green intensity values of candidate neighborhood region f 8−9 S Area , S ConA Area and convex area of candidate region f 10 S Sol Ratio of the area to the convex area f 11 S Ext Ratio of pixels in the candidate region to pixels in the total bounding box f 12 S P er Distance around the boundary of the region f 13 S CirD Diameter of a circle with the same area as the region f 14−15 S AxiA , S AxiB Major and minor axes lengths of the ellipse with a same normalized 2 nd central moments as the region f 16 S Ecc Ratio of distance between the foci and the major axis length of the ellipse with a same 2 nd moment as the region f 17 S Eul Number of objects in the region minus the number of holes in those objects  (c) ellipse that has same second moments as the region (in green) with minor and major axes (in red) and ellipse's foci are shown with blue stars. filters are suitable for the detection of convex shapes and objects with a limited range of sizes regardless of their contrast with respect to the background. The CI filters evaluate the convergence degree of gradient vectors within a local area (support region) towards a pixel of interest [31].
Given an input image I(x, y), for each pixel with spatial coordinates (x, y), the convergence index (CI) is defined by where M is the number of points in the filter support region S, and ϕ (x, y, θ i , m) is the orientation angle of the gradient vector at the polar coordinate (θ i , m) with respect to the line, with direction i, that connects (θ i , m) to (x, y). The angular difference ϕ, is given by where α is the image gradient orientation within the convergence filter support region. The support region polar coordinates are denoted by the radial coordinate m, the distance from point of interest (x, y) in pixel, and angular coordinate θ i which is sampled with N equally spaced radial lines The set of radial lines is emerging from the point where the filter is being applied to, and is equally distributed over a circular region centered at the point of interest (x, y).
Several CI filters have been proposed according to the way how the support region is defined [31]. Among different CI filters, the adaptive ring filter (ARF) [32], the sliding band filter (SBF) [33] and super-elliptical filter (SEF) [34] are more suitable for the MA detection, since they can be parameterized to use a narrow band for spotting the donut-shaped objects in the gradient-weighted images.
As shown in Fig. 7a, the ARF has a ring-shaped region of support and its radius changes adaptively. The response of the ARF is obtained via: where N is the number of support region lines as described previously, d corresponds to the width of the ring (band), and R max represents the outer limit of the band.
The result of applying equation 7 on the input image I is the filter's response image. For each candidate point we can obtain the radius of support region at that location. The shape estimation is performed by searching for the radius of the ring support region for each candidate in the image: where R ARF is the radius of the support region that corresponds to the highest convergence for location (x, y).
The SBF support region shown in Fig. 7b is a band of fixed width with varying radius in each direction, where the maximization of the convergence index at each point is obtained by: where R min and R max represent the inner and outer sliding band limits, respectively.
Given the more flexible shape formulation for the SBF support region, the shape is defined by N independent radii: where R SBF (x, y, i) has a different value for each direction and R SBF (x, y) represents the radius average. The super-elliptical band in the SEF filter allows the model to characterize a larger variety of shapes whilst at the same time reducing the irregularity in shape by its 2-fold symmetry constraints. The SEF filter is defined by r+d m=r cos(ϕ (j+i,m) ) + cos(ϕ (j−i+ N 2 ,m) ) + cos(ϕ (j+i+ N 2 ,m) ) + cos(ϕ (j−i+N,m) ) , To simplify the equation, we used ϕ (i,m) instead of ϕ(x, y, θ i , m) which represents the angle between the gradient vector at point (i, m) and the direction that is currently being analyzed. In order to consider the possible orientations of the super-elliptical filter, the parameter j is introduced in this filter as shown in Fig. 7c. The shape estimation in the SEF is more complex than in the case of the ARF and SBF filters. Here, the R SEF (x, y) is calculated by averaging the major and minor radii of the obtained super-ellipse.
The described ARF, SBF and SEF filters are applied on the normalized image (I N ) separately, and the filters' responses and estimated radii at the center of each candidate are included in the set of features. The same process is repeated by applying the filters on the gradient-weighted image (I W OS ). In total 12 LCF-based features are extracted as described in Table I. The local convergence index filters' responses and corresponding estimated shapes and radii on a normalized image patch (I N ) are shown in the first row of Fig. 8 while the second row illustrates the results on the gradient-weighted image (I W OS ).
The complete feature set F for the final classification step includes 29 features, which are 7 intensity-based features, 10 shape-based descriptors and 12 LCF-based features as shown in Table I. For each candidate in the M set, all the features are extracted and the entire feature set is rescaled to a normal distribution with zero mean and unit standard deviation.

E. Supervised Classification
To discriminate MAs from non-MA candidates, we use a hybrid sampling/boosting algorithm, called RUSBoost, proposed by Seiffert et al. [35]. RUSBoost is an adaptive boosting classifier (AdaBoost) [36] in combination with a random undersampling technique (RUS) [37] which is specifically designed to improve the performance of models trained on imbalanced data.
Boosting algorithms combine weak classifiers into a single ensemble classifier through majority voting. In particular, the AdaBoost algorithm [36] starts by setting all instances' weights equally. Then, in each iteration, a weak model is formed by the base learner and then based on the calculated error the weights are adjusted in a manner that the weights of misclassified instances are increased while the weights of correctly classified instances are decreased. Weak models are added sequentially in each iteration and trained using the weighted data to correctly classify previously misclassified instances. The process continues until a defined number of weak learners have been created or no further improvement can be made on the training set.
Random undersampling (RUS) is a technique to deal with the imbalanced class problem by changing the class distribution of the training data set [37]. The sampling is performed simply by randomly discarding instances from the majority class until a specific class distribution is obtained. In particular, the RUSBoost randomly undersamples a subset from the majority class in each iteration of the AdaBoost algorithm and consequently, the weak learners are trained using a balanced training set. For new input data, each weak learner generates a prediction value, which is weighted by the learners stage value. The final class is assigned by the summation of all weighted prediction values.
Within the context of MAs classification, RUSBoost (with decision trees as the weak learners) is a suitable classifier since we deal with a skewed set with the minority of MA candidates and the majority of non-MA candidates.
The details of RUSBoost classifier are given in Algorithm 2. The algorithm takes as input a training set (F 1 , L 1 ), ..., (F m , L m ), where F i ∈ F is a feature vector and l i ∈ L is its label. The weight for each sample i in iteration t is denoted by W t (i), where initially all weights are set to 1/m. In each iteration, a temporary class-balanced training set D is obtained by randomly undersampling the D set. At each round, the weak learner algorithm searches for h t that minimizes the error with respect to the selected subset D and distribution weights W . Then the weights are updated by calculating the overall error t and the weight update parameter α t . The new weight distribution W t+1 is used to train the next weak learner, and the process is iterated T times. The final strong classifier H(x) is a weighted combination of T weak classifiers which gives us the probability of being MA for candidate x. Fig. 9 shows examples of a RUSBoost classifier result in which the probability of being MA for each candidate is Psedo-loss calculation for D and Wt: Update weights and normalization: W t+1 (i) = Wt(i) exp(−αtl i ht(F i ))/ m i=1 W t+1 (i) end for Output and final classifier: H(x) = T t=1 αtht(x)/ T j=1 αt indicated using a heat color map. The yellow circles represent the annotation by experts.

A. Material
For the evaluation of proposed method, we used three publicly available datasets, i.e. the e-ophtha-MA [38], Di-aretDB1 [39] and Retinopathy Online Challenge training set (ROC) [40] and also our two new public RetinaCheck datasets called RC-RGB-MA and RC-SLO-MA [41], [42]. The specifications of these datasets are summarized in Table II including the imaging modality, image resolution, camera angle view and number of annotated MAs.
1) e-ophtha-MA: e-ophtha is a public database of color fundus images designed for scientific research in Diabetic Retinopathy [38]. It contains 233 healthy images (i.e. no lesion) and 148 images with microaneurysms or small hemorrhages manually annotated by ophthalmologists. The images have four different resolutions, ranging from 1440 × 960 to 2544 × 1696 pixels with 45 • field of view (FOV).
2) RC-RGB-MA: The RC-RGB-MA is a retinal image dataset which is collected in the framework of the RetinaCheck project managed by Eindhoven University of Technology, the Netherlands [41]. The 250 RGB images in this dataset are acquired with a DRS non-mydriatic fundus camera with a resolution of 2595×1944 and 45 • FOV. Two experts annotated the MAs in all images using the Microaneurysm Annotation Tool (RC-MAT) [43]. The agreement values between two experts are shown in Table III. This dataset is used to assess the potential of the proposed method as a standalone application in a large-scale DR screening setting i.e. RetinaCheck project.
3) RC-SLO-MA: In addition to the RGB fundus images, we also provide a public dataset called RC-SLO-MA in which the images are captured using the Scanning Laser Ophthalmoscopy (SLO) technique [42]. The images in this dataset are acquired with an EasyScan camera (i-Optics Inc., the Netherlands) using both green and infrared lasers. The RC-SLO-MA dataset includes 58 images with a resolution of 1024 × 1024 and a 45 • field of view. The MAs are once annotated by an expert using only the green images and another expert manually labeled the MAs using the infrared laser images. As we can see in Table III, the agreement value between two experts on the labeled MAs is less than the value obtained for the RC-RGB-MA dataset. On the other hand, the agreement between two experts on the discrimination between images with MAs and images without MAs is higher using the SLO images compared to the RGB ones. The overall agreement results show that MA detection is a challenging task for the human experts as well.
4) DiaretDB1: The DiaretDB1 is a publicly available dataset comprises 89 color fundus images [39]. Four medical experts marked the MAs independently and reported confidence levels {< 50%, 50%, 100%} which are representing  the certainty of the decision that a marked finding is correct (Table III). Because of the disagreement between the four experts annotations, a consensus of agreement higher than 75% is used to assign an MA label to a region (resulting in 182 MAs). There are 3 to 5 dark spot shaped artefacts caused by a dirty camera lens which are located in the exact same position in the images of the DiaretDB1 dataset.

5) ROC:
The Retinopathy Online Challenge (ROC) [40] contains 50 training images and 50 test images, where the annotations for the training set are publicly available. Since for the test images the gold standard is not provided and a new submission to this challenge is no longer possible, we use only the 50 training images to train and test with a 10fold cross validation approach.

B. Parameters Optimization
As we described in Section II-E, in the final classification step we use the RUSBoost algorithm with the decision trees as its weak learners. The implementation requires the determination of a number of parameters including: 1) Maximum number of splits: a parameter required for the construction of the decision trees. For a maximum number of splits equal to 1, each node represents a stump. By increasing the number of splits, the complexity of the decision trees will increase. 2) Learning rate λ: a shrinking parameter which reduces the contribution of each weak learner in order to prevent overfitting. The optimal choice of the learning rate depends critically on the number of trees; a small λ requires a large number of trees to achieve good performance [44]. 3) Number of trees: is the number of trees used for training.
Increasing the number of trees results in a decrease in both the bias and variance terms in the bias-variance decomposition. On the other hand, using too many trees may result in overfitting [44]. To investigate which parameters would be optimal, we trained classification ensembles for up to 5000 trees, varied the maximum number of splits from 1 to 512 and used as learning rates 0.1, 0.25, 0.5, 1. The optimization is done on 10% of e-ophtha-MA samples using all extracted features. A part of this experiment is shown in Fig. 10, demonstrating the mean squared error (MSE) for the maximum number of splits equal to 32, 64, 128, 256, 512.
The optimal parameters obtained in this experiment were 4727 for the number of trees, 128 as the maximum number of splits and a learning rate of 0.5. It is clear from Fig. 10 that for a large enough number of splits, the RUSBoost classifier does not seem to lead to overfitting, even in case λ = 1.00. Using too many weak learners and a high number of splits leads to a long computation time but a small decrease in the MSE. The curves in Fig. 10 show a very small decrease in the MSE values by increasing the number of trees higher than 1000. Therefore, we use 1000 as the number of trees, 100 as the maximum number of splits and 0.5 as the learning rate.

D. Microaneurysm Detection Evaluation
For the evaluation of proposed MA detection, we performed repeated 10-fold cross-validation for each dataset separately. In this approach, each dataset is divided randomly into ten equally sized partitions. Each partition is used as test data, while the other 9 partitions are used for training the classifier. The cross-validation procedure is repeated 10 times, yielding 10 performance results which are then averaged to produce a single estimation. Examples of MA detection on the five datasets are given in Fig. 11.
To measure the performance of MA detection, we used the free-response operating characteristic (FROC) curve [45] by plotting the sensitivity against the average number of false positives per image (FPI). Sensitivity represents the proportion of MAs correctly detected by the algorithm, while FPI is the number of non-MAs wrongly detected as MAs. For the sake of comparison with the other methods, the sensitivity values for the FPI rates values of 1/8, 1/4, 1/2, 1, 2, 4, and 8 were obtained from the FROC curve. The final FROC score (F score ) is defined by the average of sensitivity values at these seven predefined FPIs [40]. In addition, we also obtained the partial area under the FROC curves (F AU C ) between 1/8 and 8 FPI using trapezoidal integration and normalization by dividing with the maximum FPI [17].
The FROC curves of the proposed MA detection method on the five datasets are shown in Fig. 12a. Table V compares the F score and the F AU C of the proposed method with the state-of-the-art on the e-ophtha-MA, ROC and DiaretDB1 datasets. Fig. 12b illustrates the results on the RC-RGB-MA dataset where the classifier is trained three times using the annotations provided by expert 1, expert 2 and the agreement of both experts, respectively. The performance of the two human experts are also demonstrated in this figure.

E. Image Classification Evaluation
The proposed MA detection method extracts several candidates per image and assigns a probability value to each candidate. For the image classification, the highest probability among all candidates is considered as the image score. If the image score is greater than a certain threshold level, it will indicate that the image has signs of diabetic retinopathy and contains at least one detected MA. A low image score shows the absence of MAs indicating a healthy retina. For the valuation, the area under the ROC curve (AUC) is obtained for all five datasets as shown in Fig. 12c in which the human expert performance on the RC-RGB-MA dataset is also demonstrated by a red plus sign.

F. Features Importance Analysis
In the section, we investigate the importance of the extracted features to show the relative contribution of different types of features for the MAs classification. The feature importance analysis for the RUSBoost classifier is performed by summing changes in the risk due to splits on every predictor and dividing the sum by the number of branch nodes. This sum is taken over best splits found at each branch node. The feature importance associated with this split is computed as the difference between the risk for the parent node and the total risk for the two children. The nodes are split based on the impurity which is dependent to the split criteria. Here we calculated the importance values using the following split criterion: 1) Gini's Diversity Index [47]: is the observed fraction of classes with class i that reach the node. A pure node (with just one class) has Gini index 0; otherwise the Gini index is positive.
where p(i) is defined the same as for the Gini index, and a pure node has deviance 0; otherwise, it is positive.
3) Twoing rule [47]: where L(i) denotes the fraction of members of class i in the left child node after a split, and R(i) denotes the fraction of members of class i in the right child node after a split. P (L) and P (R) are the fractions of observations that split to the left and right respectively. If the expression is large, the split made each child node purer. Similarly, if the expression is small, the split made each child node more similar to each other, and hence more similar to the parent node, and so the split did not increase node purity. We evaluated all the 29 features using the Gini's Diversity Index, Deviance and Twoing rule on a subset of images from the e-ophtha-MA dataset. The obtained feature importance maps are shown in Fig. 13 and the corresponding description of each feature is given in Table I. The features with the top 12 maximum importance values in Fig. 13 are selected as the reduced subsets, which are then used to train the RUSBoost classifier. The intensity-based features and the proposed local convergence filter-based features are frequently used in decision making for all the three split criteria, while the shapebased features have the least contributions.
We also trained the classifier using each category of features individually or combined with each other. The results are shown in Table VI     Although the intensity-based features have the most substantial contributions (F score = 0.366 and F AU C = 0.473) compared to two other feature categories (see Fig. 14), inclusion of the LCF-based features improves the performance significantly, resulting in a F score of 0.518 and a F AU C of 0.604. On the other hand, including shape-based features improves the performance only slightly. As illustrated in Fig. 13, the most important descriptor in the category of intensity-based features is the minimum green intensity (G min ), and among the LCFbased features the response of the SEF filter on the gradientweighted image (F W SEF ) has the highest contribution.  [20], [24] using WebPlotDigitizer application [46].

IV. DISCUSSION AND CONCLUSION
The proposed method is in the form of a pipeline of techniques where the performance of each step depends on the output of the previous step. Hence, the method is validated extensively in each one of three stages: candidate extraction, MA detection, and image classification.

A. Candidate Extraction Performance
The goal of the candidate extraction is to reduce the computational burden by decreasing the number of objects for further analysis in the next step. As shown in Table IV, the proposed candidate extractor outperforms the state-of-the-art methods by achieving a sensitivity value of 0.82 on the ROC dataset. The average of false positives per image is higher than the ones reported by other methods, however, it is less than 0.06% of the total number of pixels in one image. In this step, it is important to include true MAs as much as possible without considering the balance between sensitivity and specificity. The false positives are later discarded in the classification step. The evaluation on other datasets demonstrates that the introduced candidate extractor performs better on the e-ophtha-MA (Sen = 0.95) and the RC-RGB-MA (Sen = 0.94) since both datasets contain high resolution and high quality images.

B. Microaneurysm Detection Performance
The detection of MAs highly depends on the imaging device characteristics and image properties such as resolution, modality, compression technique, illumination and contrast variation. For this reason, we evaluated the proposed method on five datasets. The results are compared with the state-ofthe-arts on the e-ophtha-MA, ROC and DiaretDB1 datasets, while the RC-RGB-MA and RC-SLO-MA datasets are used for the comparison with human experts and for the validation on a different image modality (SLO images).
On the high quality images of e-ophtha-MA dataset our method achieves a F score of 0.510 and a F AU C of 0.575 which are significantly higher than the values reported by Wu et al. [23]. On the same dataset Zhang [19] used contextual descriptors for the classification of MAs with F score equal to 0.440 and with a slightly better partial AUC of 0.586. Fig. 12b and Table V show the performance of the proposed method on RC-RGB-MA, as well as the inter-expert variability between two human experts. The F score values of 0.534 and 0.614 are obtained by using the annotations provided by the 1 st expert and 2 nd expert respectively while using the agreement of two experts on the annotations results in a better performance (F score = 0.647). In all three cases of different training approaches, the proposed method has a similar or slightly better performance than the human experts.
The RC-SLO-MA dataset provides a unique possibility for the validation of the proposed method on a different image modality (SLO). The green laser usually provides a better contrast between vessels and background, while the infrared laser penetrates deeper into the retina and visualizes pathologies in different layers of the retina. The proposed method performs slightly better on the green channel than the infrared image while the obtained F score values for both channels are relatively lower than the values achieved on the RGB datasets since the SLO images typically contain more background noise than conventional fundus images (RGB) and have a lower spatial resolution.
The images of the ROC dataset were acquired at different resolutions using different cameras. The variability in image resolution as well as the presence of noise and artifacts makes it more challenging to detect MAs in the images of this dataset. However, our method overcomes these difficulties and achieves a F score of 0.471 outperforming the state-of-theart approaches. The sensitivity values at predefined FPIs are compared with the others in Table V showing our method has higher sensitivities at FPI values of 1/8, 1/4 and 1/2. Wang et al. [25] achieved slightly higher sensitivity at 1 FPI. Although some methods [17], [18], [22], [24], [48] reported higher sensitivities at 2, 4 and 8 FPI, it should be noted that the FPI of 1.08 is considered as an indication of "clinically acceptable" FPI and higher values are not adequate for a computer-aided diagnosis system in clinical environments [40].
Although most of the images in DiaretDB1 dataset contain several dark spots and artefacts, our proposed method achieves remarkably better results compared to other methods. As shown in Table V, our method achieves a F score value of 0.547 and a F AU C of 0.564 for the DiaretDB1 dataset. These values are significantly higher than the ones based on the dynamic shape features presented by Seoud et al. [24] and the gradient vector analysis technique proposed by Dai et al. [22]. Fig. 12a illustrates the ROC curves of image classification, where the highest probability value among all candidates is assigned to the image as its score. The highest AUC is obtained for the e-ophtha-MA dataset (AUC = 0.96) while the lowest one is for DiaretDB1 dataset (AUC = 0.89). The dark spot artefacts in the DiaretDB1 have high probability values which can cause a wrong image classification resulting in a lower AUC. Although the MAs detection on the RC-SLO-MA dataset has a lower performance than the other datasets, the high AUC value of 0.91 is obtained for the image classification. Visual inspection, as well as the high AUC, reveals that the obvious MAs are much easier to detect in SLO images compared to RGB images. The comparison between the ROC curve of our method on the RC-RGB-MA dataset (AUC = 0.90) and the sensitivity/specificity of the 2 nd expert (reference: 1 st expert ) demonstrates the relatively high success rate of the proposed method in DR detection (see Fig. 12a).

D. Computation Time
The fully automated MAs detection method is developed in MATLAB 2016b (MathWorks, Inc.) with an average computation time of 3 minutes per image using an Intel Core i7-5820 CPU at 3.30 GHz.

E. Conclusion
In this paper, we have proposed a new method for detecting MAs in retinal images using a gradient weighting technique, a new set of features based on local convergence filters (LCF) and a random undersamping boosting classifier. Feature importance analysis demonstrates that LCF-based descriptors can well characterize the low contrast MAs since the LCF filters are based on gradient convergence and not intensity. The performance of the proposed method on the ROC, DiaretDB1 and e-ophtha-MA datasets shows the competitiveness of the introduced approach against state-of-the-art techniques. Moreover, the evaluation results on five public datasets demonstrate that the proposed MAs detection method is insensitive to the characteristics of the imaging device, image resolution and image modality. Future work will involve exploiting the introduced LCS-based features for the detection of dot hemorrhages and bright lesions namely exudates.