A combined fuzzy backtracking search optimization algorithm to localize retinal blood vessels for diabetic retinopathy

Fuzzy c-means Clustering (FCM) algorithm allows to categorize different pixels in an image by employing fuzzy memberships [27]. Let P = {P₁, P₂, P₃,...,P_M } represents M pixels in an image. These M pixels are to be segregated into C clusters by minimizing the fitness function $({ \mathcal J })$

$$\begin{eqnarray}&&{ \mathcal J }=\displaystyle \sum _{q=1}^{M}\displaystyle \sum _{r=1}^{C}{u}_{{rq}}^{\mu }{\parallel {P}_{q}-{V}_{r}\parallel }^{2};C\leqslant M\end{eqnarray} \tag{ 1 }$$

where u_rq depicts fuzzy membership value of qth pixel in rth cluster. V_r and μ represents rth cluster center and fuzziness control parameter respectively.

The membership value describes the likelihood value of a particular pixel to a cluster center. The likelihood is purely computed based on the distance between pixel and the respective cluster center. The membership value for a particular pixel is high if the pixel is close to cluster center and vice-versa.

The cluster center (V_r ) and membership function (u_rq ) are upgraded by using equations (2) and (3) respectively.

$$\begin{eqnarray}&&{V}_{r}=\displaystyle \frac{{\sum }_{q=1}^{M}{u}_{{rq}}^{\mu }{P}_{q}}{{\sum }_{q=1}^{M}{u}_{{rq}}^{\mu }}\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&{u}_{{rq}}={\left[\displaystyle \sum _{s=1}^{C}{\left(\displaystyle \frac{\parallel {P}_{q}-{V}_{r}\parallel }{\parallel {P}_{q}-{V}_{s}\parallel }\right)}^{2/(\mu -1)}\right]}^{-1}\end{eqnarray} \tag{ 3 }$$

Initially, the cluster centers are chosen randomly and are iterated until values of similar pixels belongs to a particular cluster [27].

Backtracking search optimization algorithm (BSA) is an adaptive evolutionary algorithm developed for determining global minimum or maximum value [71]. BSA consists of both historic and evolution populations. Historic population stores the randomly chosen members in earlier generations for creating a search direction matrix. While the evolution population contains the present information of the members. The following steps briefly describes BSA.

2.2.1. Initialization

Let us assume that (N, D) be the size of population and dimension of the variable, respectively. Then, the evolution (H) and historical population (H^old ) are initialized randomly as given in equation (4) and 5, respectively.

$$\begin{eqnarray}H=\left[\begin{array}{ccccc}{H}_{11} & {H}_{12} & {H}_{13} & \ldots & {H}_{1D}\\ {H}_{21} & {H}_{22} & {H}_{23} & \ldots & {H}_{2D}\\ \vdots & \vdots & \vdots & \ldots & \vdots \\ {H}_{N1} & {H}_{N2} & {H}_{N3} & \ldots & {H}_{{ND}}\end{array}\right]\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}{H}^{{old}}=\left[\begin{array}{ccccc}{H}_{11}^{{old}} & {H}_{12}^{{old}} & {H}_{13}^{{old}} & \ldots & {H}_{1D}^{{old}}\\ {H}_{21}^{{old}} & {H}_{22}^{{old}} & {H}_{23}^{{old}} & \ldots & {H}_{2D}^{{old}}\\ \vdots & \vdots & \vdots & \ldots & \vdots \\ {H}_{N1}^{{old}} & {H}_{N2}^{{old}} & {H}_{N3}^{{old}} & \ldots & {H}_{{ND}}^{{old}}\end{array}\right]\end{eqnarray} \tag{ 5 }$$

Values in H and H^old are distributed uniformly and lies in the interval of $\left[{H}_{\min },{H}_{\max }\right]$.

2.2.2. Selection-I

In this step, a new historic population is generated for evaluating the search directions. The rules for generating the new historic population is given as

$$\begin{eqnarray}{H}_{{ij}}^{{old}}=\left\{\begin{array}{lc}{H}_{{ij}} & {{ifa}}_{1}\lt {a}_{2};i=1,2,..,N;\forall j=1,2,..,D\\ {H}_{{ij}}^{{old}} & {otherwise}\end{array}\right.\end{eqnarray} \tag{ 6 }$$

where a₁ and a₂ are random numbers generated between interval $\left[0,1\right]$. The rows of H^old matrix are randomly shuffled by using a permutation function.

$$\begin{eqnarray}&&{H}^{{old}}={permuting}({H}^{{old}})\end{eqnarray} \tag{ 7 }$$

2.2.3. Mutation

A trail population (T) is generated by computing the difference between historic (H^old ) and evolution population (H). The equation for generating the T is

$$\begin{eqnarray}&&T=H+{ \mathcal F }\left({H}^{{old}}-H\right)\end{eqnarray} \tag{ 8 }$$

where, ${ \mathcal F }$ is the scale factor, which regulates the mutation step length. The historic knowledge of the population passes through complete evolutionary procedure.

2.2.4. Crossover

In this step, the final trail population (T^final ) is computed. T^final is evaluated by using equation (9).

$$\begin{eqnarray}{T}_{{ij}}^{{final}}=\left\{\begin{array}{ccc}{H}_{{ij}} & {if}\quad {{map}}_{{ij}}=1 & ;i=1,2,3,\ldots ,N\\ {T}_{{ij}} & {otherwise} & ;j=1,2,3,\ldots ,D\end{array}\right.\end{eqnarray} \tag{ 9 }$$

where, map is a binary integer valued matrix and leads the crossover direction. A detailed discussion to generate map matrix is available in [71]. The values of T^final are regenerated again if they fall outside the boundary limits.

2.2.5. Selection-II

A set of new individual population for subsequent iterations are generated in this step. Greedy selection algorithm is employed to select the individual population. The new individual population is computed as

$$\begin{eqnarray}{H}_{i}=\left\{\begin{array}{ccc}{T}_{i}^{{final}} & {if}\quad { \mathcal J }({T}_{i}^{{final}})\leqslant { \mathcal J }({H}_{i}) & ;i=1,2,3,\ldots ,N\\ {H}_{i} & {otherwise} & ;\forall j=1,2,\ldots ,D\end{array}\right.\end{eqnarray} \tag{ 10 }$$

where ${ \mathcal J }(.)$ is the fitness evaluation function.

To note, these steps are iterated until a global minimum or global maximum is attained [71].

Preprocessing step is performed on the green channel of the image. The green channel is chosen as it corresponds to the highest luminance and therefore the possibility of identifying large contrast between background and retinal blood vessel in this band. Further, brightness correction is also carried out to compare the global reflectance, on all over the image. This step simply determines the global mean brightness (G^mean ) over the entire image and then passes the image through a window of sufficient size to guarantee that the mean brightness inside the window $({W}_{{ij}}^{{mean}})$ is proportionate to the global mean G^mean . The updated brightness of the center pixel of the window is described as

$$\begin{eqnarray}&&{b}_{{ij}}^{{new}}=\displaystyle \frac{{b}_{{ij}}^{{old}}}{{W}_{{ij}}^{{mean}}}\times {G}^{{mean}}\end{eqnarray} \tag{ 11 }$$

where ${b}_{{ij}}^{{new}}$ and ${b}_{{ij}}^{{old}}$ are new and old pixel values.

The objective function $({ \mathcal J })$ of fuzzy clustering is minimized by using BSA. The hybrid algorithm is used to avoid local minima problem. Fuzzy-BSA segmentation is illustrated in below steps.

• a.  
  The green channel corrected brightness image is considered for clustering. The initial parameters of FCM such as cluster number (C = 3) and fuzzy control parameters (μ = 2) are defined. Similarly, the initial parameters of BSA like population (N = 15), scale factor $({ \mathcal F }=0.3)$, stopping criterion and mixed rate are also described. The number of iterations of Fuzzy-BSA is considered as 75.
• b.  
  The evolution matrix (H) is randomly initialized as described in equation (4).
• c.  
  For individual cluster center in (H), fitness values is evaluated by using equations (1), (2) and 3 respectively.
• d.  
  Begin BSA iteration and in each iteration calculate the fitness function $({ \mathcal J })$ value and revise the fitness and H.
• e.  
  Stop BSA if stopping criterion is encountered and export global minimum value and final u matrix.
• f.  
  Create C clustered images by utilizing final u_rq matrix.

In the post-processing stage, small, individually connected components are removed, small gaps are filled, and connected components with low thinness ratios are joined. With size 3 × 3 and rank 7, a rank order filter was used to fill in the slight gaps and remove the small connected components. By using a thinness ratio lower than a predetermined threshold, the binary connected components are eliminated[58]. The closure operator with size 2 × 2 is used to fix the small holes and border distortion. The equation for thinness ratio (TR) is given as

$$\begin{eqnarray}&&{TR}=\displaystyle \frac{4\ast \pi \ast \mathrm{Total}\ \mathrm{area}\ \mathrm{of}\ \mathrm{individual}\ \mathrm{connected}\ \mathrm{components}}{\mathrm{Total}\ \mathrm{perimeter}\ \mathrm{of}\ \mathrm{individual}\ \mathrm{connected}\ \mathrm{components}}\end{eqnarray} \tag{ 12 }$$

Structured Analysis of Retina (STARE) [38] and Digital Retinal Images for Vessel Extraction (DRIVE) [12] are used to evaluate the performance of proposed Fuzzy-BSA technique. To note, DRIVE database comprises of 40 color fundus images, which are saved in JPEG compressed format. Out of 40 retinal images, 7 images contains exudates, pigment epithelium changes and hemorrhages. The images are acquired in Netherlands from a diabetic retinopathy screening camp. Samples of abnormal and normal images from DRIVE database is depicted in figure 2. And as can be seen, the uncertainty between retinal vessels and the background particularly for the abnormal images can be clearly visualized. In addition to that, the white and dark spots produce heterogeneous brightness over the entire image.

Zoom In Zoom Out Reset image size

STARE database, on the other hand, consists of 20 eye fundus color images in which 10 images comprises pathology. TopCom TRV-50 fundus camera was used to capture images in STARE database with an FOV of 35°. The resolution of each image is 700 × 605 pixels with 8 bits of quantization. All the images are available in Portable Pixmap Format (PPM) format. To note, two sets of ground truth image segmentation are available which are manually done by two separate ophthalmologists. However, only the second manual segmentation is widely opted by the researchers for evaluating the performance. Samples of normal and abnormal images from STARE database is shown in figure 3. It is evident from this image that the complications are similar to that in DRIVE database. However, more number of abnormal images are present in STARE dataset.

Zoom In Zoom Out Reset image size

In addition to this, we also collected a dataset containing 2500 retinal images from various hospitals in Bangalore and Mangalore. Of them, the major contributors are Nethradhama Super Speciality Eye hospital and A J Institute of Medical Sciences and Research Center. The dataset contains both normal and abnormal retinal images. The ground truth images are drawn by Dr. Eshwar Fani from AJ Institute of Medical Sciences with the help of Dr. Ramachandra Shet.

The proposed algorithm is executed in Intel Core i5 CPU with speed 2.30 GHz and 8 GB RAM. The code is written in MATLAB 2021b version. Specificity, sensitivity and accuracy are the primary performance metrics considered for evaluating the proposed algorithm. To note, these measures are computed for individual images from each dataset and final average values are compared against the state-of-art algorithms.

$$\begin{eqnarray}&&{Specificity}=\displaystyle \frac{{TN}}{{FP}+{TN}}\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{Sensitivity}=\displaystyle \frac{{TP}}{{FN}+{TP}}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{Accuracy}=\displaystyle \frac{{TN}+{TP}}{{FN}+{TP}+{TN}+{FP}}\end{eqnarray} \tag{ 15 }$$

where, TP True Positive (TP) represents that the segmented pixels belongs to ground truth. True Negative (TN) denotes that the non-vessel pixels correctly segmented as non-vessel pixels. False Negative (FN) symbolizes the vessel pixels segmented as non-vessel pixels, whereas False Positive (FP) represents non-vessel pixels segmented as vessel pixels. Specificity determines the probability of the segmentation algorithm that will correctly recognize non-vessel pixels. Similarly, sensitivity determines the probability of the segmentation algorithm that will correctly recognize retinal blood vessel pixels. Finally, accuracy renders the overall performance of the segmentation algorithm.

Figures 4 and 5 depicts the results of pre-processing stage for DRIVE and STARE dataset. Figure 2(a) and (b) shows the normal and abnormal retinal images of DRIVE dataset respectively. Similarly, figure 3(a) and (b) illustrates the normal and abnormal retinal images of STARE dataset respectively. Figures 4(a), (d), 5(a) and (d) shows the green channel images of the original images. Similarly, figures 4(b), (e), 5(b) and (e) depicts the illumination correction of figures 4(a), 4(d), 5(a) and (d) respectively. To note, preprocessing step is homogeneous for both the normal and abnormal datasets which therefore rectifies the brightness of the image while maintaining contrast between background and blood-vessels (see for instance figures 4(c), (e), 5(c) and (e)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 6(a), (d), 7(a) and (d) represents the output of Fuzzy- Backtracking Search optimization algorithm (Fuzzy-BSA). In general, the fuzzy clustering algorithm tries to minimize the intra cluster spectral variations therefore tiny (diameter) gray pixels are omitted as the general clustering algorithm fails to differentiate between gray color and the surrounding pixels. But, the backtracking search optimization algorithm improves the segmentation result by concentrating even on tiny gray pixels. The BSA targets to maximize the thinness parameter which results in attaining the cluster centers (towards the thin connected units) with homogeneous spectral values for the blood vessels even if it infringes the inter cluster variable measure. The BSA optimization has a distinct impact on the final image depending on whether or not the tiny blood vessels shows originally in Fuzzy cluster.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 6(b), (e), 7(b) and (f) depict the post processing results for the output of Fuzzy-BSA stage. Each individual pixels are removed by using rank-order filter. Distortions in boundary and tiny holes are rectified by using closing operation. However, in figure, 7(f) the area of pigment epithelium is more and creates a ambiguous content with the retinal blood vessels. In turn, this makes the proposed algorithm a challenge for segmentation and fails at this condition.

Figures 8(a), (c) and figure 8(b), (d) shows the normal and abnormal images collected from the hospitals and with corresponding segmented outputs, respectively.

Zoom In Zoom Out Reset image size

Table 2 summarizes the performance metrics of the proposed algorithm and also compares with different state-of-art algorithms for both DRIVE and STARE datasets. '−' indicates that the authors of the corresponding papers have not computed the metrics. It is evident from this table that the proposed algorithm shows superior performance when compared to the state-of-art algorithms with sensitivity of 0.782, specificity of 0.985 and accuracy of 0.973 for DRIVE dataset. Similarly, performance metrics of STARE dataset are sensitivity of 0.84, specificity of 0.972 and accuracy of 0.964. It can be observed from table 2 that the results of proposed algorithm are comparable with the supervised as well as unsupervised algorithms.

Table 3 depicts the performance measure of the proposed Fuzzy-BSA algorithm with the first observer values for the medical dataset collected from various hospitals. We note that the first observer is a doctor itself from A J Institute of Medical Sciences and Research. The performance metrics of Sensitivity, Specificity and Accuracy are 0.681, 0.852 and 0.894 respectively.

Figures 9(a)–(c) shows the comparison of proposed Fuzzy-BSA algorithm with second observer values for DRIVE, STARE and medical datasets.

Zoom In Zoom Out Reset image size

The robustness of proposed method can be observed for both the normal and abnormal images. The results are identical even with images having hemorrhages, exudates and change in pigment epithelium for DRIVE database. Similarly, merely constant results are obtained for STARE database which has more number of abnormal images. No varying results are found even while performing Monte-Carlo runs due to applying of two level optimization technique.

Table 4 depicts the comparison of proposed method execution time with different retinal vessel segmentation techniques. The proposed algorithm needs an average of 52s to segment an retinal blood vessel image and thus making it compatible as compared to widely used techniques in the literature. We note, the execution time can be further reduced by implementing batch processing.