Algorithmic Analysis of Vesselness and Blobness for Detecting Retinopathies Based on Fractional Gaussian Filters

: All around the world, partial or total blindness has become a direct consequence of diabetes and hypertension. Visual disorders related to these diseases require automatic and specialized methods to detect early malformations, artifacts, or irregular structures for helping specialists in the diagnosis. This study presents an innovative methodology for detecting and evaluating retinopathies, particularly microaneurysm and hemorrhages. The method is based on a multidirectional Fractional-Order Gaussian Filters tuned by the Differential Evolution algorithm. The contrast of the microaneurysms and hemorrhages, regarding the background, is improved substantially. After that, these structures are extracted using the Kittler thresholding method under additional considerations. Then, candidate lesions are detected by removing the blood vessels and fovea pixels in the resulting image. Finally, candidate lesions are classiﬁed according to its size, shape, and intensity properties via Support Vector Machines with a radial basis function kernel. The proposed method is evaluated by using the publicly available database MESSIDOR for detecting microaneurysms. The numerical results are summarized by the averaged binary metrics of accuracy, sensitivity, and speciﬁcity giving the performance values of 0.9995, 0.7820 and 0.9998, respectively. H.A.-R. and J.G.A.-C.; Data curation, J.M.C.-D., I.C.-A. and J.R.-P.; Visualization, M.d.J.E.-A., I.C.-A. and J.M.C.-D.; Writing—original draft preparation, H.A.-R., and M.d.J.E.-A.; Writing—review and editing, J.G.A.-C. and J.M.C.-D. and J.R.-P.; Funding acquisition, J.G.A.-C. and J.R.-P. All authors have read agreed the published version


Introduction
The stress, inappropriate feeding behaviors, lack of interest in prevention methods, besides the extensive periods without proper medical screening, have increased the spread of retinal maladies in the world. About 126 million people suffered a kind of Diabetic Retinopathy (DR) in 2011. Such a number could increase up to 191 million by the end of 2030, with around 56.3 million people suffering a vision-threatening complication [1,2]. Regular medical screenings on vulnerable patients (e.g., persons suffering from diabetes or hypertension) may detect ocular problems to prevent permanent blindness in up to 98% of the cases. Although these programs have resulted in a declining prevalence and incidence distinctive shape, besides the proximity to the optic nerve, fovea, or blood vessels make difficult the accurate recognition of the suspicious lesions.
This work aims aims to contribute to the accurate recognition of retinopathies, particularly microaneurysm and hemorrhages. For that purpose, we present an innovative framework powered by a multidirectional filter based on the Fractional-Order Gaussian Function (FOGF). This filter is adjusted via a well-known optimization algorithm, such as Differential Evolution (DE). The application of such a filter meaningfully improves blood vessel contrast, as well as positively enhances the delimiting region of retinal aneurysms. The framework also includes the Kittler thresholding method under additional conditions for extracting the microaneurysms and hemorrhages structures. Thenceforth, candidate lesions are detected by removing the blood vessels and fovea pixels from the resulting image. Finally, these lesions are classified according to their size, shape, and intensity properties via SVM with a radial basis function kernel. The proposed method is evaluated by using the publicly available database MESSIDOR for detecting microaneurysms. Results show astonishing characteristics of this method for identifying microaneurysm structures in most of the studied cases.
The remainder of this document is organized as follows: Section 2 briefly describes the medical concepts related to retinopathy diagnosis. Section 3 presents the proposed methodology, as well as some practical considerations. Subsequently, Section 4 details the methodology carried out and the corresponding experimental results. Finally, Section 5 summarizes the most relevant conclusions.

Background
Retinography is a noninvasive diagnostic technique that certainly does not use contrast agents in the acquisition process [17]. This technique produces color images from the inner part of an eye, well-known as Fundus Images (FIs), which can be used to detect some diseases, e.g., diabetic retinopathies and glaucoma. These images are characterized by a quasi-circular shape with a small rectangular flange that represents the Field of View (FoV) of the capture device. Figure 1 shows the typical internal structures composing a fundus image.

Optic Disc Macula
Fovea Blood vessels In general, the fundus images from healthy-patients have four common elements: the optic disk, blood vessels, fovea, and macula (see Figure 1). The Optic Disk (OD) is a circular and bright structure that interconnects the internal and external parts of the eye. Crossing the OD, there are the blood vessels and optic nerve [8]. The blood vessels are typically reddish tubular structures that exchange oxygen and move essential nutrients from and into the human eye. These structures, identified as veins and arteries, are spread in the image as a complex network connected with smaller elements known as venules and arterioles in their extremities. Found close to the optic disk is the fovea, visually recognized as an oval or circular red region. This element is extensively irrigated by the choroid, which is a membrane free of blood vessels. Located in the fovea, the biological photoreceptor cells, namely cones and rods, transform the light stimuli into electrical signals transmitted by the optic nerve to the brain. Additionally, the fovea is surrounded by the macula that is specialized in the perception of details. Under certain circumstances, specific features of these elements are used to determine the physical or mental health of a patient. For instance, recent studies have shown a direct link between the retinal state and its association and potential use in the diagnosis of dementia [18], Parkinson's [19], and Alzheimer's diseases [20], as well as other cognitive deficits [21] and neurodegenerative conditions [22]. In the same context, modifications in shape and optic disk size are a clear pointer to serious illness like glaucoma. However, different circumstances influence the appearance of extraneous biological structures in FI, most of them attributed to retinopathies. Retinopathy is a noninflammatory retina disease that can also include other medical conditions [10,18,20]. Particularly, diabetic and hypertensive retinopathies are the most studied and well-known retinal diseases. Both medical conditions may provoke similar symptoms (i.e., extraneous structures in the image), but the number, shape, size, and apparition probability may distinctively characterize each disease. Besides, its presence in the retina is highly influenced by the disease evolution degree. These representative structures are used as a starting point for the diagnostic assessment of retinopathies. Figure 2 illustrates the most distinctive structures characterizing retinopathies, such as neovascularization, hemorrhages, exudates, and microaneurysms.  Specifically, the neovascularization (NV) is related to a new blood vessel generation produced by a prolonged lack of blood in specific regions of the eye [23]. It is usually a consequence of diabetic retinopathy. Nevertheless, due to the stressing process of how these new blood vessels were created, they generally have irregular, weak, and low-quality walls prone to rupture. Plus, these blood vessels may come up as saliencies in the contour of the optic nerve, as shown in Figure 2a. The retinal hemorrhages (H) typically appear when blood vessel membranes get weak and start to bleed. But it is also provoked by changes in blood composition or disturbances in circulation [24]. Hence, blood is accumulated in the retina or vitreous humor. These hemorrhages produce irregular reddish spots in the fundus image, like those observed in Figure 2b. In critical cases, hemorrhages may produce partial or complete blindness because photoreceptor cells suffer from light occlusion and are barely stimulated. Moreover, the exudates (E) are typically generated when proteinic substances are spread into the eye from blood vessels, forming brilliant depots varying their size from small spots until large zones in the FI. Some exudates surrounding the fovea are displayed in Figure 2c. Like NV, this pathology is also caused by diabetic retinopathy. The last characteristic structure of retinopathy is the microaneurysm (mA), which appears as small reddish spots in different zones of the FI (see Figure 2d). These could be easily confused with small punctual hemorrhages. However, the mA is induced by the blood vessel dilations but not by draining blood out to the vitreous body. This specific set of undesired structures is used for diagnosing both the Retinopathy Degree (RD) and Macular Edema Risk (MER). Furthermore, Table 1 contains the respective constraints for the diagnostic evaluation of RD and MER, considering the occurrence frequency of microaneurysms, hemorrhages, and neovascularization. Table 1. Specific constraints to evaluate the Retinopathy Degree (RD) and Macular Edema Risk (MER), by employing the ocurrence frequency of microaneurysms (mA), hemorrhages (H), and neovascularization (NV), as well as the MaCula to Exudate distance ( MCE min ) and the optic disk diameter (D OD ) [25].

Level
Retinopathy Degree (RD) Macular Edema Risk (MER) Under the diagnosing parameters, a patient is considered free of any retinopathy when RD = 0, while a patient suffering the worst symptoms of such disease may get RD = 3. Moreover, MER is evaluated considering the occurrence of exudates; when they are detected, there is a substantial risk of macular edema if the MaCula to Exudate distance is inferior to the optic disk diameter, MCE min ≤ D OD [26]. Consequently, the number of mA and H can be used, almost independently, as an efficient estimator of retinopathies. Such a premise is the central axis of this study.

Materials and Methods
The proposed framework uses images from the public database MESSIDOR [25]. This database comprises 1200 images aimed to develop methods for diagnosing diabetic retinopathies. Each image has the RD and MER diagnosis provided by experts according to the disease evolution. For evaluation and learning purposes, this work employs the mA delineations generated by Habib et al. in [27,28]. The flowchart depicted in Figure 3 presents the main blocks of the proposed method.  At first, the input color image is decomposed into the corresponding RGB channels. The green channel is selected because it presents the highest contrast, instead of computing the traditional gray level images by averaging the three RGB channels. This choice is computationally more efficient, and according to color theory, the green channel contains more details and less noise than other channels. In the upper-branch of the processing flow, the Fractional Order Gaussian Filters (FOGF) and Kittler thresholding method are applied to the selected image channel. The FOGF filters are proposed to ease the task of detecting the small H and mA structures. These operate mainly as feature descriptors and functionally enhance the detection of the structures by improving the contrast. Thus, the filter parameters are tuned by using a metaheuristic optimization method. The Differential Evolution (DE) algorithm was chosen due to its efficiency and relative simplicity of implementation. Hence, the FOGF filters are focused on detecting H and mA. Moreover, given their fractional-order nature, they are applied in eight preferential orientations. Next, a thresholding process based on spatial constraints over the Kittler minimal error method is applied to the filtered image, which extracts several complex structures from the background. In the lower-branch of the process, an algorithm for identifying both blood vessels [29] and fovea based on an efficient conics detection method [30] was implemented. A differential outcome is obtained by subtracting the partial results from both branches (see Figure 3) to discriminate the structures of interest. At this stage, some pixels belonging to other structures such as the fovea or blood vessels are found. These are eliminated by identifying the pixels related to a specific structure and then removing them from the binarized image. As a result, a set of candidate lesions is found. Afterward, a series of five features are extracted from each candidate lesion, i.e., the size, eccentricity, mean, and minimum intensity of the FOGF filter response, and length of the major axis of the ellipse that contains the candidate lesion. Finally, from these features, mA and H lesions are found by applying an SVM classifier over this collection of features.

Fractional Order Gaussian Filters
Since the successful application in anomalous diffusion [31], the fractional calculus has found a vast number of uses on pattern recognition, image denoising [32,33], and texture enhancement [34].
In the literature, there exists a variety of definitions for the fractional differential operator, such as Liouville, Riemann-Liouville [35], Caputo, Caputo-Fabrizio [36], Grünwald-Letnikov [37], Weyl, Marchaud [38], Hadamard [39], Chen [40], Chen-Marchaud [41], Fourier transform [42], among others [43]. In practice, this operator is selected according to its properties and suitability to the problem to analyze. So, we implemented the operator presented by Tseng et al. [44], which is defined as follows is the Fourier Transform of the same signal. This definition is used to generate the FOGF filter bank.
Thus, four filters are shown in Figure 4 by using the Gaussian function N x (µ, σ) as the base function f (x), since µ is the mean and σ 2 is the data variance. Using Fractional-Order Gaussian Filters as a texture descriptor is pretty recent in image processing. Moreover, these are turning out especially suitable for texture applications when fractional orders ν ≤ 1 are employed. In fact, a recent study of the fractional calculus applied in image processing was extensively discussed in [45]. Such work presented the discretization of most used fractional-order derivatives validating their advantages on improving image analysis. Hu et al. proposed a fractional differential operator mask with a non-integral adaptive step and fractional order to enhance the texture and analyze its features [34]. In practice, it is noticeable that FOGF presents similar behavior as the Difference of Gaussians (DoG) filters, which may extend the number of applications.

Model Parameters Selection
All the FOGFs are generated using the kernel D ν x N x (µ, σ) obtained with (1) and the normal distribution function N x (µ, σ). In such a kernel, a set of four parameters should be estimated, i.e., the fractional-order ν, the domain x, the mean µ, and the standard deviation σ. However, to determine x two parameters (n and χ) are required, such that x is rendered by n elements K i uniformly distributed in the interval [−χ, χ]. Figure 5 shows an illustrative example of the kernel parameters. It is noteworthy that asymmetry is detected on this profile, which is produced by the fractional nature of the filter. Withal, the vector of kernel parameters is represented by X = (ν, χ, n, µ f , σ f ) having six design elements, i.e., X ∈ R 6 . Determining the optimal values of these six parameters requires an exhaustive evaluation that is generally time-demanding and impractical. In particular, if each parameter is tuned at a time. Naturally, to obtain the reliable and optimal performance of the proposed filters, it is fundamental to simultaneously determine all the parameters and anticipating a certain degree of dependency between them. Therefore, the vector X is estimated by using the Differential Evolution (DE) algorithm.

Differential Evolution (DE) Algorithm
DE is a well-known optimization method characterized by its relative tuning simplicity and effectiveness for solving multiparametric problems. This algorithm was proposed by Storn and Price [46] as a method based on an evolutive mechanism. DE starts generating an initial population at random, and finds the optimal solution by fitting the experimental data to a given model representing the problem. This technique uses some strategies to combine the initial population for the posterior choice of the best individuals per iteration. After each combination, a mutation function is applied to generate new diversified individuals. The canonical implementation of this algorithm comprises four stages: initialization, mutation, crossover, and selection. Hence, there are N candidates in the population P, each candidate has a position (solution) vector given by X ∈ R D , which are initially estimated via the uniform distribution into the search space. During each generation (iteration) G, every individual in the population is considered as an objective vector X G i . In the mutation stage, at least two donor elements X G r i , randomly chosen with r i ∼ U I (1, N) and i {r i } = ∅, are combined to create the new donor solution vector u G+1 I that may interact with the next generation G + 1. In this work, the combination strategy is defined by where F 1 and F 2 are the mutation factors, and X G best is the objective vector corresponding to the best model adjustment in the current generation. Experimental experiences allowed choosing the values for these factors such as F 1 = 0.5 and F 2 = 0.01. For the crossover stage, some internal components of the donor vectors are randomly modified by using a threshold C R ∈ [0.0, 1.0]. The crossover function produces a test vector v G i that diversifies the population at each iteration. The components of v G i are obtained as shown where is a random number generated from the standard uniform distribution, and I u ∼ U I (1, D) is an integer randomly selected from the interval [1, D]. Moreover, the threshold value C R = 0.5 was determined from an exhaustive experimentation process. Finally, during the selection phase, the best individuals satisfying the model are chosen among the test v G i and objective vectors X G i . The last three stages of the DE algorithm are successively repeated until the convergence criterion is satisfied.

Objective Function
The performance of the response image Im R is analyzed to select the candidate vectors containing the optimal parameters. This image is obtained by adding all the responses of the retinal image Im convolved with the FOGFs bank D ν x N x (µ, σ). Given that fractional filters are intrinsically asymmetric, they were applied on eight directional configurations, i.e., θ k = 2πk/8, ∀ k ∈ {0, 1, . . . , 7} to avoid preferential directions (spatial bias). These configurations can be visually represented by {→, , ↑, , ←, , ↓, }. Such a generalization to 2D fractional filters was previously studied by Chen et al. in the context of image-enhancing [47]. Therefore, Im R is determined through where * is the convolution operator, and D ν x N k x (µ, σ) corresponds to the k-th FOGF using the parameters X G i oriented with the directional configuration k. Furthermore, it is required to define an objective function to evaluate the performance of the FOGF for a given candidate vector X G i during the DE process. Such a function is stated as since | · | is the cardinality operator.

Kittler Thresholding Method
This histogram-based thresholding method was proposed by Kittler and Illingworth [48]. It searches a minimal probabilistic decision error in a supposing bimodal distribution. They considered that a gray level image with g = 2 n pixel values can be separated into two planes, background and foreground, by using the image histogram h(g). So, h(g) represents both planes by the addition of at least two probability distributions p i (g). These functions are modeled as normal distributions using only the means µ i and variances σ 2 i . Besides, by employing the a priori probability P i ∀ i ∈ {1, 2}, the complete distribution becomes According to the Bayes theorem, the optimal threshold T h separating both distributions yields the minimal decision error. This threshold can be found by using the following reduced function where P i (T) = ∑ b g=a h(g), µ i (T) = ∑ b g=a gh(g)/P i (T), and σ 2 i (T) = ∑ b g=a (g − µ i (T)) 2 h(g)/P i (T) for the interval [a, b] established as [0, T] for i = 1 and [T + 1, n] for i = 2. Therefore, the minimal error threshold T h is estimated such as However, for images having a high density of low-intensity pixels, the obtained threshold could be useless as they are polarized to dark gray levels. Notice that the error function proposed by Kittler presents a high probability of finding points with a low-intensity when J(T|T → 0) J(T|T → ∞). Consequently, reliable thresholds are found to some extent far from the low bins in the histogram. Extensive experimentation allowed to accomplish such condition by finding the inflection point between the interval [T h , max(J(T|T > T h ))], which is depicted in Figure 6.
This rule is applied even if this point was computed as the point T i such as J (T i ) ≈ 0 and verified a sign change in J (T i ). In this work, for the sake of the numerical stability, the threshold is approximated to the averaged point T m k in the interval T ∈ [T h , max(J(T|T > T h )] as given

Blood Vessels and Fovea Correction
Blood vessel extraction and fovea detection are not the central subjects of this study, but considering their intensity features, they could appear next to aneurysms and hemorrhages after thresholding. Therefore, these structures should be identified and removed previous an extensive image analysis. In such a sense, some methodologies are then proposed for this detection.

Fovea Detection
Concerning the ellipsoidal (or circular) shape of the fovea, this structure may be correctly identified by applying a conic detector filter on the retinal image Im R . However, some lesions may appear close to the fovea leading to inaccurate results and misdetection of microaneurysms and hemorrhages. Thence, an intensity-based segmentation that preserves the fovea irregularities is used. In that context, the estimated threshold image f th is obtained by binarizing the green channel of Im R via the Kittler thresholding method. Next, the undesired tubular and punctual structures are removed by using the morphological operator having a circular shape and radius 12 (D S 12 ) structural element. The foveal image is then computed by It is essential to notice that given the problem nature, it may be necessary to use a complete version of the fovea, instead of the one given by the image opening. This element can be obtained via the region growing given by Im f as the seed over image f th . Figure 7 depicts an example of the proposed methodology to extract the fovea region.

Blood Vessels Detection
Blood vessel extraction uses the Frangi method [49]. This method employs a specialized filter to enhance tubular structures, which is based on calculating the eigenvalues λ i of the Hessian matrix Nonetheless, to control the width in the tubular structures, the direct partial derivatives are substituted by the smoothed partial derivatives of a Gaussian filter N XY = 1 √ 2πσ exp − x 2 +y 2 2σ 2 with mean µ = 0 and variance σ 2 , such as Thus, the vascularity function V o in (16) is determined by the eigenvalues of H N XY , where the variables α, β, and γ represent the thresholds controlling the filter sensibility.
Different diameter structures are detected using a variety of scales in the Gaussian filter. Hence, the result is obtained by computing the maximum value of the vascularization function sweeping over different scales. Considering a certain number N of scales in a searching range σ ∈ [σ min , σ max ], the vascularity image is calculated by For this study, the particular parameters for the vessel detection are σ ∈ [2,3] and N = 21 (i.e., steps of 0.05), and (α, β, γ) = (1.0, 2.25, 15) . Due to the nature of this image (vasculature with both high and low intensities), the conditions for the application of modified Kittler thresholding are not achieved. Therefore, the original Kittler algorithm should be used. Nevertheless, the original Kittler method also fails to find the complete vascularity network. This negative effect is seen in Figure 8b, where only the major vasculature is completely detected. Given that the main objective of this work is the removal of blood vessels in the original image rather than a precise detection, a length filter with a fixed threshold is exploited. This last threshold is applied to the normalized version of the vascularity image (V l = V l → [0, 255]). Under these circumstances, a threshold value of t c h = 0.01 is used to create the candidate V c s image by V c s =V l > t c h . Finally, the vascularity image is segmented utilizing the length filter given as where t a h corresponds to the minimal area that is maintained by the n-th region in V c s = n V c s (n) to be considered as a blood vessel. A representative result using the fixed length filter (t a h = 100) is also presented in Figure 8d.

Candidate Lesion Classification
Even though the removal of fovea and blood vessel pixels could produce a set of candidate lesions, such collection could still include undesired elements such as unconnected vessels, spurious lesions, or speckle noise. All candidate lesions are passed through an SVM classifier to decide whether a candidate could represent a microaneurysm (mA) or hemorrhage (H). It uses an optimization method to identify support vectors s i , weights w i , and bias b i , which are employed to classify vectors x into a class c according to where k is the kernel function. For this study a radial basis function kernel, k(s i , x) = exp − s i − x 2 , was used. The candidate lesions represented by a vector x is composed of five features: size, eccentricity, the average and minimum intensity of the FOGF response, and length of the major axis of the ellipse containing the candidate lesion. Each feature is selected to cope with the specific properties of undesired candidates.

Numerical Results
All experiments in this study were carried out in a PC with Intel R Core TM i5-6200U CPU @ 2.3-2.4 GHz, and 8 GB RAM. The algorithms were coded and evaluated in MATLAB R R2016a running on Microsoft R Windows TM 10 OS. The design of optimal filters for improving contrast as well as enhancing features are complex tasks, but fundamental in medical imaging. Such a complexity relies on a high number of parameters for tuning and a close correlation among them. The adverse effects of this correlation are reduced by simultaneously adjusting all filter parameters. For this purpose, the DE optimization algorithm was implemented, and its parameters C R , F 1 , and F 2 were analyzed by combining exhaustively of their values. The chosen values that provide the best results were C R = 0.5, F 1 = 0.5, and F 2 = 0.01. Figure 9 depicts the objective function value (error) evolution during the tuning process. Bear in mind that the evolution of f obj (cf. (6)) is plotted for 100 generations. An averaged error of 0.1888 is found only in the first generation. However, the DE method evolves towards a better solution with each new iteration. Therefore, the model coefficients can be selected from a wide range of combinations (or search space) without significant variation in errors. In the implemented version, the parameters were constrained to positive values and rounded to the upper integer. DE allowed determining the best performance by tuning the FOGF, as Figure 10b shows. Subsequently, Figure 10c presents the resulting image after applying these filters to a retinal image. Still, FOGFs detect not only microaneurysms but also blood vessels, fovea, and hemorrhages. Such undesired misclassifications can be attributed to the objective function properties and the intensity similarities between the target structures. Contrarily to the resulting image presented in Figure 10b, the exudates are distinguished and extracted from the one shown in Figure 10d. This last result was computed using the Kittler thresholding, where microaneurysms, hemorrhages, blood vessels, and some fovea pixels are accurately detected. It is noteworthy that the method selectivity on detecting microaneurysms and hemorrhages is boosted by using peculiar features, shapes, and sizes. Figure 11 exhibits microaneurysms and fovea detection. Previously, the fovea and blood vessels removal function was efficiently executed.  Notwithstanding, a severe retinopathy (G R = 3) was correctly diagnosed in this image. The redundancy of information and the sufficient number of the structural elements allowed obtaining this complicated diagnosing. So, four binary metrics, accuracy (Acc), balanced accuracy (BAcc), sensitivity (Sen), and specificity (Spe) were selected for evaluating the classification of proposed method [50]. These measurements are determined from a binary image A and a reference image B (ground-truth) as follows where TP stands for True Positives, i.e., points designated as positive in binary image A that also appear as positives in reference image B (or simply A B); TN corresponds to True Negatives, i.e., points in A marked as negative that also appear as negatives in B (orĀ B ); FN means False Negatives, i.e., points established as negatives in A but which appear as positives in reference image B (orĀ B); FP defines the False Positives, i.e., points labeled as positives in A that are designated as negatives in B (or A B ). B andB symbolizes the Positive and Negative counted values of B respectively, and | · | represents the cardinality operator. These metrics produce a continuous output in the interval [0.0, 1.0] since 1.0 indicates the perfect matching between the images A and B. The evaluation of microaneurysms detection is summarized in Table 2. There, an average balanced accuracy of 0.8909 and an average sensitivity of 0.7820 were singularly representatives. Note that all detected regions were considered to build this table. Accordingly, the results contain only a few pixels belonging to hemorrhages or blood vessels. Thus the size could allow discriminating between microaneurysms and hemorrhages, the retinopathy degree is assigned by computing the number of detected elements from each one of the classes. Results of microaneurysms and hemorrhages detection are shown in Figure 12. The proposed method gives the results depicted in Figure 12b, and the result after applying the Kittler thresholding method is presented in Figure 12c. In that image, the reference microaneurysms are displayed in green and those detected by the proposed method in blue. As expected, the result obtained by applying the thresholding method contains the majority of microaneurysms. Nevertheless, the proximity of these structures with some blood vessels gave rise to their elimination from the final result. In consequence, an efficient and robust method to detect and remove blood vessels is highly necessary for improving results. Moreover, most dispersed points in the reference image belong also to pixels from blood vessels, but they correspond to narrower blood vessels than the regular, continuous structures.

Conclusions
In this paper, an innovative framework for enhancing and detecting microaneurysms in fundus images is presented. It applies a multidirectional filter based on Fractional-Order Gaussian Filters (FOGFs). On the one hand, the internal parameters of this filter were simultaneously adjusted by using the well-known Differential Evolution (DE) algorithm. On the other hand, a modification to the Kittler thresholding method to cope with histograms containing a high density of low-intensity points was also implemented. The proposed method was evaluated on the publicly available database MESSIDOR for detecting microaneurysms. This methodology provided acceptable results for enhancing and detecting microaneurysms, which an average accuracy of 0.9995 and an average sensitivity of 0.7820 for the global detection process. Moreover, the method was able to enhance and detect hemorrhages that were separated from microaneurysms using shape and size features. Microaneurysms rendered irregular shapes, contrarily to hemorrhages that are characterized by punctual or quasi-circular shapes. Under these circumstances, the obtained results can be utilized to classify the degree of retinopathy present in patients by employing as support the parameters given in Table 1. The number of detected structures is highly relevant for microaneurysms and hemorrhages diagnosing. Thus, when a single microaneurysm point is detected, this should be included in the general account directly affecting the diagnosis. Finally, the proposed method is extremely competitive for detecting at least one microaneurysm structure in most studied cases. Finally, the proposed approach can found some potential applications for detecting similar patterns in medical imaging. For instance, on digital mammograms, the effective localization of microcalcifications is highly demanded, in such a problem, the early detection of tiny abnormal structures in the woman's breasts is a priority. Moreover, the microcalcifications are also commonly found in affections derived from medical disorders in the brain, liver, and thyroid gland, which justify a broader prospective evaluation of the proposed method.

Conflicts of Interest:
The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations and symbols are used in this manuscript: