CFLHCF: Simultaneous Detection of the Optic Disc and Exudates Using Color Features, Local Homogeneity and Contextual Features

In this section, the proposed CFLHCF and MMBS‒SED for EXs detection is developed, which following five stages: 1) image pre-processing, 2) elimination of the OD, 3) extraction of EXs via color component features, 4) selected samples of EXs, and 5) EXs classification, respectively. The fundus images are modified in the first stage to normalize illumination, contrast enhancement, noise removal, color space selection, and optic disc localization. In the second stage, CFLHCF is adopted to obtain EXs candidates. Next, MMBS‒SED is used to distinguish EXs from non‒EXs candidates. The following subsections discuss the details of each stage.

2.1 Image preprocessing

We apply a histogram specification technique [20] to the Red (R), Green (G), and Blue (B) chanel images to the original fundus images. The intensity of the original images is normalized and equalized with the reference histogram images selected by human experts. Segment reallocation and intersegment transformation are used to obtain the transformed histograms for both images with and without edges. These two histograms utilized to perform histogram specifications are combined to get the desired histogram. In this example, the input image’s gray level, i, is mapped to another gray level, d, by Eq. (1):

$C_{\text {in }}(i)=C_{\text {desired }}(d)$                      (1)

where, C_in(i) and C_desired(d) are the cumulative distribution and histogram, respectively, derived from the original fundus image. In other words, we look for the gray level, d, by Eq. (2):

$d=C_{\text {desired }}^{-1}\left(C_{i n}(i)\right)$                        (2)

We use Eq. (2) to transform each gray level in the original image into a normalized image. The output image is transformed back to RGB if the input image was originally RGB. However, at this stage, as observed from the image’s center out to its periphery, the brightness gradually reduces when the distance to the center image is increased (see Figure 1(b)). To increase the contrast between the EXs lesions in the image, a local contrast enhancement approach [21] must be used. The local contrast method was applied as follows: Given a small N×N running window w and the intensity of a color component image for the pixel range p, then the color image is converted to the full range image f_new [0, L‒1] with the linear stretching using Eq. (3):

$f_{{new }}=(L-1)\left(\frac{\psi_w(p)-\psi_w({Min})}{\psi_w({Max})-\psi_w(\operatorname{Min})}\right)$                        (3)

where, ψ(Min) and ψ(Max) are the sigmoidal function of the minimum and maximum value of the enhanced features among RBG color component of the fundus image p, and L=256 for 24‒bit images. The sigmoid function, which is also known as the standard logistic function is generated by applying Eq. (4):

$\psi_w(p)=\left[1+c x p \frac{\left(m_w-p\right)}{s_w}\right]^{-1}$                          (4)

In Eq. (4) derivatives are given as Eq. (5) and Eq. (6) and applied to all images to contrast windowing:

$m_w=\frac{1}{N^2} \sum_{(i, j) \in w(k, l)} \quad p(i, j)$                       (5)

$s_w=\sqrt{\frac{1}{N^2} \sum_{(i, j) \in w(k, l)}\quad \left(p(i, j)-m_w\right)^2}$                         (6)

where, m_w and S_w are the input image’s local mean and standard deviation-correspondingly, which changes for every location. As usual, w(k, l) and N are the locations of the pixel within the window and the number of window sizes. Considering Eq. (4), the gain value determines the real contrast by controlling the speed at which the function changes from minimal to maximum. For the upper limit of acceptable standard deviation values and window size, we select S_w=0.2 and N=59×59, respectively. While the local enhancement method increases the luminosity of images, it also increases the brightness of noise or artifacts, which can lead to incorrectly classifying these pixels as EXs lesions. Hence, the third pre-processing step used image denoising by a median filtering method [22]. The image was filtered with an arithmetic mean filter of size 3×3, the smallest filter required in this instance to remove the most visible signs of impulse noise. Finally, we adopted the J model selection criterion to measure various color model efficiency to select the suitable model color selection criterion [23]. The matrix function J calculated the within‒class (S_w) and between‒class (S_b) scatter matrices. The class reparability of the region of interest pixel class for various color models, Eq. (7), can be used to get this function.

$J=\operatorname{trace}\left(\frac{S_b}{S_w}\right)$                        (7)

The parameter J is a set of overall clustering color pixels, while the denominator is assumed to be the number of regions where an image should be segmented. Here, S_w is the distribution of sample points around each mean vector, and S_b is the mixing weights of the sample pixels around the mean vector. The method for estimating both clusters can be defined as Eq. (8) and (9) respectively.

$S_b=\sum_{i=1}^C N_i\left(M_i-M\right)^2$ with $M=\frac{1}{N} \sum_{i-1}^C N_i M_i$                          (8)

where, C is the estimated color pixels regarded as one class and N_i is the image that has N pixels at the ith color pixel class C_i. The whole samples in the image form dataset are $N=\sum_i N_i$.

$S_w=\sum_{i=1}^C S_i$ with $M_i=\frac{1}{N_i} \sum_{n \in C_i}^C\left(X_n-M_i\right)^2$ and $M_i=\frac{1}{N_i} \sum_{n \in C_i} X_n$                              (9)

where, C is the number of clusters, X_n is the sample color points n, and M_i is the mean vector of class C_i. Based on the suitable color model selection, six delicately selected models are used in this experimental testing: RGB, YIQ, LUV, LAB, HSL, and HSI.

In the experiments, we used 89 fundus images to test the effect of color model selection on the function J selection criterion. The quantitative results for the different color models are shown in Table 1. A higher value of J indicates quantitative results in a more detailed segregation classification, while other members within each class are closer. Based on the experimental testing, the final choice of color model is the LUV color model. Therefore, this stage converts the RGB color space into LUV and evaluates the processed luminance space as a combination with empirically set weights.

Table 1. A comparative analysis of color models

2.2 Elimination of the optic disc (OD)

Due to the candidate EXs’ similarity to a yellow region in fundus images in terms of shape, color, intensity, brightness, and sharpness, the elimination of OD is a critical step in the proposed methods to classify the EXs [24, 25]. The area in the fundus image with the brightest yellow lesion is the OD. At other times, the OD is not discernible in the yellowish component and has to be localized in the red component. Therefore, OD elimination is not an easy matter. Some of these difficulties in OD elimination are shown in Figure 2.

In this work, the method attempted to locate OD accurately by a mathematical operator and threshold technique. In order to remove blood vessels from the OD region and approximate the OD, we first used the elimination method utilizing mathematical based on the closing operator (Figure 3(b)). Object contours are smoothed, thin connections are broken, and thin protrusions in the OD region are removed while morphologically closing is used on EXs. As shown in Figure 3(b), the results of closing with a structuring element weight of 7.8 were used. Afterward, the resultant image was binarized by automatic thresholding according to the OTSU threshold [26], see Figure 3(c). The result of the thresholding was a binary image where the OD found region was marked with 1’s and the background of the fundus image with 0’s. To take advantage of this feature, all the OD pixels from Figure 3(c) were used to create the candidate OD region in Figure 3(d). The largest pixel did not correspond to the true OD pixels, which is the step’s fundamental weakness. As a result, a mathematical dilation operator with a flat disk shape structuring element weight of 7.5 was applied to extend the OD region. The dilation achieved with this element is shown in Figure 3(e) (the OD region is larger to accommodate the dilation). Then, the OD was binarized using thresholding with a weight of 0.65 (Figure 3(f)). The associated OD components were included in the new binary image, as was previously mentioned. The “good” of these results was finally localized. The OD could then be classified as the fundus image’s largest circular connected component. It was discovered that the OD region was cropped as ∼3,200 pixels in the fundus image. With the help of the OD, elimination could process only the pixels of the OD and omit the background pixels by applying Eq. (10) [24, 25].

$C=4 \Pi \frac{ { Area }}{ { Perimeter }{ }^2}$                        (10)

where, C is the ratio of the OD region’s area to its perimeter square’s compactness, Area denotes the size of the OD candidate area, and perimeter represents the number of pixels along its perimeter (see Figure 3(g)). Then, a binary mask of the boundary candidates is overlaid on the original image (see Figure 3(h)).

2.jpg

Figure 2. Structure of OD and varying size, color, intensity, and location in fundus images

3.jpg

Figure 3. Output images of the OD elimination. (a) pre-processed image, (b) applying morphological closing operator, (c) thresholded image using OTSU threshold method, (d) OD is eliminated on the pre-processed image, (c) applying morphological dilation operator, (f) thresholded image using OTSU threshold method, (g) candidate OD regions using Eq. (10), (h) final OD elimination and superimposed on the original image

2.3 Extraction of exudates via color features and local homogeneity

When a fundus image is provided, function J from the preceding section converts the original image into LUV space. At this point, the fundus image was expected to contain P pixels, and C initialized cluster centers. The metric interval (M) is calculated as $M=\sqrt{P / C}$.  In order to calculate the distances between pixels and cluster centers, color information values in the potential minor neighborhood should be in a 3×3 representation of the EXs region. After then, each cluster center is updated using the identified pixel centroid. This procedure is repeated until there is no longer any change in the distance between successive cluster centers. In Color Features and Local Homogeneity (CFLH), the distance (D) between the cluster centers and the color R, G, B or L, U, V pixel positions is calculated [27]. The cluster center of pixel distance from pixel x to y is defined by Eq. (11) and (12).

Distance $=\sqrt{D_{ RGB }^2+\left(\frac{ m }{ M }\right)^2} D _{ XY }^2$                         (11)

Distance $=\sqrt{ D _{ LUV }^2+\left(\frac{ m }{ M }\right)^2} D _{ XY }^2$                        (12)

where, the weight factor, m, has to have a range of 0 to 255 (RGB) and 0 to 1 (LUV), M is the region size, and x and y are the color values of the EXs pixel at the position of x and y in the fundus image. In CFLH, several regions were extracted as features, including blood vessels, red spots, microaneurysms, and related to non‒EXs. The EXs extraction is based on color information with in-depth features to discriminate between EXs and non‒EXs components. The feature of EXs lesions is that they both appear as yellow dots and are more visible than the surrounding tissue. The extraction was necessary for collecting and providing data on EXs lesions in fundus images because the authors and supervisors didn’t have any medical expertise in fundus images when the research began. Three professional ophthalmologists were able to highlight EXs locations on fundus images effectively. The EXs regions should be marked (see Figure 4); we can select the class of the EXs.

4.jpg

Figure 4. Marking EXs regions with an expert ophthalmologist

The whole EXs regions are represented by yellowish areas in the color component features. Hence, these regions and color shading were targeted for classifying the EXs regions. After being careful with the color model selection in the previous section, we compared RGB and LUV color information that were well‒suited for classifying the whole EXs‒region. The color information values in the possible minor neighborhood are the size of 3×3 used to represent the EXs region in each model, as shown in Table 2.

Table 2. The color information values to represent the EXs‒region in the smallest possible neighborhood is of size 3×3

In order to distinguish EXs candidates from background feature candidates, a set of color information was used to prove that the yellowish color has relativity with the EXs‒region in terms of the intensity feature. In this way, the EXs pixels in RGB and LUV channels (see Table 2) that the EXs have different non-EXs contrast in other channels.

2.4 Candidates of exudates using contextual features

In the previous section, the above color component features analysis uses yellow to represent the region of the whole EXs while red, brown, and orange colors represent non‒EXs regions. An example of the proposed novel contextual feature is illustrated in Figure 5. Consider Figure 5, yellow pixels are the position of the current EXs candidate region in the RGB channel, and the corresponding neighbor EXs positions are marked as blue circles.

5.jpg

Figure 5. Example of a manually cropped image containing EXs regions of interest, (a) original fundus images; (b) color information, and a close‒up view of the EXs‒region

Essentially, the proposed contextual feature is an adjustment to the mean RGB value of each candidate. Therefore, in order to distinguish EXs (yellow features) from other color features, the RGB distance between neighboring pixels is taken into consideration. In RGB color information, the probability color values range from 0 to 255 (R:0‒255, G:0‒255, B:0‒255). Therefore, the minimum (Min) and maximum (Max) intensity of the pixels values within the RGB (red, green, and blue) color model are 0 (no color) and 255 (full color). The LUV color model, in contrast, is used to describe the color’s brightness or intensity from 0 to 1. The function is chosen as the local mean and global mean value inside each color information, creating the non‒EXs and EXs regions that are used in the proposed novel contextual feature to calculate the average color value in RGB and LUV color images. Here, let L_mean represent the local mean intensities of the color channel, Con_feature denotes a novel contextual feature defined as Eq. (13) and (14), respectively.

$L_{\text {mean }}=\frac{\text { Total intensitics value in the window }}{\text { Number of window size }}$                        (13)

Con $_{\text {Feature }}=\frac{\text { Total value of the local mean }}{\text { Total of window }}$                       (14)

To understand that a local mean intensity creates even ranges between the Min and the Max values (the weighted average values), Table 2 shows the relationship between the intensities value in the window and the number of window sizes. Consider in red (R) color channel in an image, (226, 225, 221, 224, 224, 220, 220, 219, 216), the sum of that point’s nine neighbors divided by 9. To obtain the weighted average of maximum values, the total value of 1,995 divided by 9 is 221. This weight of 221 is considered the weighted average of the EXs in the Red (R) color channel, with the smallest possible neighborhood being the size of 3×3. The neighborhood’s center is then transferred to the next nearby location. The process is repeated to produce the subsequent value of the output image in each window. In this case, the weight of 221 depends on the square of size 3×3 at the red color channel. The neighborhood’s average pixel intensity is transformed into an intensity (also called local homogeneity). Similarly, the total value of the local mean intensity with the same metric can be defined as a novel contextual feature by Eq. (12). For example, if the total value of the weighted average is 221, 221, 212, 198, and 201 in the whole window of 5. A novel contextual feature by the total value of 1,053 divided by 5 is 210.

This technique, sometimes called global means intensities. Therefore, this single value represents a novel contextual feature for classifying the EXs regions in the red space. If the confidences with a novel contextual feature of RGB color channels are 210, 189, and 201 while LUV color channels are 0.90, 0.46, and 0.99, respectively. In terms of the classification of EXs as depicted in Figure 6. Examining effects, as Figure 6(a) depicts, a novel contextual feature gets good classification results of RGB with “imag016.png” from the DiaretDB0, and Figure 6(b) the segmentation results in LUV space.

2.5 Exudates classification using MMBS-SED

Following the extraction of the candidate EXs region in the previous phase, the classification errors with various thresholds are calculated. Thresholds ranging from 0 to 1 are listed at this stage. The threshold with minimum classification error is adopted. As for the EXs, the color feature, local homogeneity, and contextual features can be regarded as the main landmark for distinguishing EXs from the other features in fundus images, which plays an important role in retinal image analyses. The proposed CFLHCF and MMBS‒SED methods mainly rely on two assumptions. First, locate the EXs appear as a yellow and bright region using optimal threshold and mathematical morphology; secondly, the EXs shape exhibited using Sobel edge detector. The EXs detection method consists of the following two parts:

1) Optimal threshold and mathematical morphology to locate the EXs regions.

2) Candidate EXs regions selection using Sobel edge detector.

Step 1: Locate the EXs by computing the histogram probabilities of the fundus image by Eq. (15) [26].

$P(i)=\frac{ { number }\{(r, c) \mid {image}(r, c)=i\}}{(R, C)}$                        (15)

where, R is the number of rows and C is the number of columns, and P_(i) is the histogram probabilities. R represents the index for a row and C is a column of the image.

Step 2: Calculate within‒class variance of EXs and non‒EXs pixels. The within‒class variance defines by Eq. (16(a)-(b)).

$\sigma_{ {Within }}^2=w_B(t)^* \sigma_B^2(t)+w_E(t)^* \sigma_E^2$                      (16a)

and

$w_B(t)=\sum_{i=0}^{T-1} p(i)$                       (16b)

where, $\sigma_{ {Within }}^2$ is the weighted sum of the variances of each cluster, $\sigma_B^2(t)$ is the variance of the pixels in the below threshold (Background), $w_{E}(t)$ is the variance of the pixels in the above threshold (EXs), $w_{B}(t)$ is the weight of background pixels, and $w_{E}(t)$ is the weight of EXs pixels respectively.

Step 3: Compute between‒class variance of EXs and non‒EXs pixels by the optimal OTSU method. This method is far faster than the simple threshold. The between-class variance is defined as Eq. (17).

$\begin{aligned} \sigma_{{ Between }}^2 & =\sigma^2-\sigma_{W { ithin }}^2 \\ & =w_B(t)^*\left(\mu_B(t)-\mu\right)^2+w_E(t)^*\left(\mu_E(t)-\mu\right)^2 \\ & =w_B(t)^* w_E(t)^*\left(\mu_B(t)-\mu_E(t)\right)^2\end{aligned}$                           (17)

where, $\sigma^2$ and μ are the overall mean and combined variance of the image, respectively. An example of the optimal threshold value of 0.68 is selected, and the result is shown in Figure 6(c). However, selecting a higher threshold value than 0.68 does not give a better result since dark EXs regions are now considered to belong to the non‒EXs.

Step 4: Mathematical morphology to locate the EXs regions. As seen in the example in Figure 6(c), the threshold method does not give adequate segmentation that provides good results. Thus, a mathematical morphology approach based on the dilation operator should be applied in this stage. Mathematical techniques are divided into larger regions based on predefined criteria. Additional criteria increase the power of the segmentation method, such as region, size, and likeness between the EXs pixel and non‒EXs. The mathematical morphology and the size of the structure element are determined as follows:

In the first round, all-white pixels cannot be marked as EXs regions or their neighborhood pixels. Therefore, in the second round, the white pixel is dilated as the flat structuring element’s size increases. In the third round, the white pixels are dilated by the size of the flat structuring element is 7.4, so dilation is performed. In the fourth round, a mathematical erosion is used to remove small pixels in the dilated mask. In this round, the size of the flat structuring element in mathematical erosion is the same as that used in dilating. In the fifth round, to classify all regions of the EXs in the mask, the algorithm computes the area of the EXs in the sixth round. In the seventh round, if the compactness is higher than a present threshold, the pixel is considered EXs regions in the eighth round. Otherwise, the pixel is regarded as a part of the non‒EXs region in the ninth and tenth rounds. If all EXs pixels have been marked, the processing task is stopped. After the EXs regions have been found, binary masks are given as input to the next stage of boundary detection.

Step 5: The boundaries of the EXs are detected at this stage using four different edge detectors, including Robert, Sobel, Prewitt, and Canny edge detectors. The Sobel operator [28] uses a 3×3 neighborhood based on a gradient operator, the convolution masked by two kernels, and the Robert edge detector [28] calculates the gradient operator around the central pixel using a 3×3 neighborhood. The Prewitt edge detector [28], and the second derivative zero-crossing point correspond to the Canny edge detector [29], with a negligible risk of weak edges. The experiment suggests a Sobel edge detector is best for detecting EXs boundaries (see Figure 6(d)). Compared to the other three edge detectors, it can also identify true EXs boundaries with the minimum number of errors.

6.jpg

Figure 6. CFLHCF and MMBS‒SE detection results of imag016.png from DiaretDB0. (a) the detection results in RGB space, (b) the detection results in LUV, (c) a thresholding result with 0.68, (d) applying Sobel edge detector

In this subsection, the experiments based on three freely accessible datasets are used to evaluate the proposed method. In these experiments, the method is assessed using pixel-based evaluation criteria. In pixel-based criteria, the connected component pixel validation is used to count the number of correctly identified pixels. Based on the above descriptions, the SEN, SPEC, and ACC values were calculated by pixel-based evaluation on, DiaretDB0, DiaretDB1, and STARE datasets under varying regions and sizes. Below is a description of the details.

4.1 The performance of the proposed method

The technique was developed using MATLAB 2019b on a desktop with CPU Intel(R) Core(TM) i7‒6700K, 4.00 GHz, and 8 GB RAM. The performance of the proposed method has achieved overall datasets with RGB and LUV color models in Table 3. The DiaretDB0 offers an SEN rate of 98.12%, a SPEC rate of 98.08%, and an ACC rate of 98.10%, respectively. In order to test the suggested method on the input dataset, we chose 297 fundus images. The experimental results on the DiaretDB1 dataset demonstrate that EXs classification SEN, SPEC, and ACC were 98.41%, 98.29%, and 98.32%, respectively. Finally, we used the proposed methods to test the STARE for EXs detection. The detection of EXs was performed using the proposed method and yielded SEN, SPEC, and ACC of 98.79%, 99.81%, and 99.34%, respectively. The proposed methods are excellent for EXs classification with the LUV color model in STARE because they obtained the best sensitivity and specificity values. However, STARE’s fundus images are smaller than those in the DiaretDB0 and DiaretDB1 databases. The excellent segmentation results in fundus images “img013.png” and “img016.png” from the DiaretDB0 and “img016.png” from the DiaretDB1 with RGB and LUV color models, which includes large and small EXs regions as shown in Figure 7. An example segmentation results in fundus images “img0306.png” and “img0359.png” from the STARE with RGB and LUV color models, which include large and small EXs regions shown in Figure 8.

7.jpg

Figure 7. Visual examples of the proposed method for EXs detection on DiaretDB0 and DiaretDB1 datasets, where the regions in black line represent EXs, (a)‒(a4) EXs segmentation results of fundus image “img013.png” from the DiaretDB0 dataset with RGB, (b)‒(a4) EXs segmentation results of fundus image “img013.png” from the DiaretDB0 dataset with LUV, (c)‒(c4) EXs segmentation results of fundus image “img016.png” from the DiaretDB1 dataset with RGB, (d)-(d4) EXs segmentation results of fundus image “img013.png” from the DiaretDB1 dataset with LUV

4.2 Comparison with the state-of-the-art method

An efficient method to segment and detect EXs in fundus images is proposed. There are algorithms in the literature for segmenting EXs using different techniques. Pereira et al. [36] offered double thresholding and Ant Colony algorithms for classifying EXs in fundus images. The proposed method was relatively good in results on HEI‒MED: the SEN of 80.82% and SPEC of 99.16%, Naqvi et al. [37] have proposed scale-invariant feature transform and K‒means clustering for classifying EXs on DiaretDB1. They reported an SEN of 92.70% and a SPEC of 81.02% for classifying EXs. Annunziata et al. [38] used a green color model for detecting EXs. They reported 71.28% SEN and 98.36% on SPEC for EXs. Imani and Pourreza [39] have proposed dynamic thresholding and mathematical morphology method for classifying EXs in fundus images on the DiaretDB1. They reported 89.01% on SEN and 99.93% on SPEC. Kaur and Mittal [40] have used iterative clustering in their algorithm to find EXs on the STARE. They reported 96.41% SEN and 96.57% SPEC with 81 fundus images. Mo et al. [41] have combined the probability map and maximal probability region for identifying true hard EXs. They reported 92.55% on SEN and 71% on predictive capabilities, but the system was tested with only 169 fundus images. The comparative analysis of the detection of EXs with the state‒of‒the‒art methods was comprehensively summarized, as shown in Table 4.

Table 3. Comparison of exudates detection performance on three different publicly available datasets

Databases	Performance measures (%) in RGB space			Performance measures (%) in LUV space
	SEN	SPEC	ACC	SEN	SPEC	ACC
DiaretDB0 (130 image)	89.01	88.90	88.04	98.12	98.08	98.10
DiaretDB1 (89 images)	90.89	89.94	89.16	98.41	98.29	98.32
STARE (78 images)	91.80	91.23	91.40	98.79	99.81	99.34

8.jpg

Figure 8. Visual examples of the proposed method for EXs detection on STARE datasets, where the regions in black line represent EXs, (a)‒(a4) EXs segmentation results of fundus image “img0306.png” from the DiaretDB0 dataset with RGB, (b)‒(a4) EXs segmentation results of fundus image “img0306.png” from the STARE dataset with LUV, (c)‒(c4) EXs segmentation results of fundus image “img0359.png” from the DiaretDB0 dataset with RGB, (d)‒(d4) EXs segmentation results of fundus image “img0359.png” from the STARE dataset with LUV

Table 4. Comparison of the obtained EXs detection performance of state‒of‒the‒art methods with the proposed method

Image processing methods have been proposed to diagnose DR in fundus images. In the future, we will combine different lesion segmentation methods based on machine learning optimizations and double thresholding [42] to classify DR into cotton wool spots, microaneurysms, hemorrhages, and EXs. If the final detection result of this combination is sufficiently good, it is possible to automate the early screening of EXs in fundus images. The automated screening application would reduce experts’ workload since only fundus images, identified as normal or abnormal by the automatic application, need to be further examined by experts.