Multi-parametric optic disc segmentation using superpixel based feature classification

Accepted Manuscript

Multi-parametric Optic Disc Segmentation Using Superpixel Based Feature Classification Zaka Ur Rehman, Syed S. Naqvi, Tariq M. Khan, Muhammad Arsalan, Muhammad A. Khan, M.A. Khalil PII: DOI: Reference:

S0957-4174(18)30770-X https://doi.org/10.1016/j.eswa.2018.12.008 ESWA 12349

To appear in:

Expert Systems With Applications

Received date: Revised date: Accepted date:

26 April 2018 3 December 2018 4 December 2018

Please cite this article as: Zaka Ur Rehman, Syed S. Naqvi, Tariq M. Khan, Muhammad Arsalan, Muhammad A. Khan, M.A. Khalil, Multi-parametric Optic Disc Segmentation Using Superpixel Based Feature Classification, Expert Systems With Applications (2018), doi: https://doi.org/10.1016/j.eswa.2018.12.008

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • A new preprocessing pipeline is used for vessel removal and optic disc enhancement • Image regions are characterized by statistical and textural properties • A set of highly discriminative features is adapted

AC

CE

PT

ED

M

AN US

CR IP T

• A comparison of supervised classifiers for optic disk localization is presented

1

ACCEPTED MANUSCRIPT

Multi-parametric Optic Disc Segmentation Using Superpixel Based Feature Classification Zaka Ur Rehman Department of Computer science and IT, The University of Lahore, Gujrat campus,Gujrat, Pakistan [email protected], T:+92 305 689 4791

Syed S. Naqvi

Tariq M. Khan

CR IP T

Department of Electrical and Computer Engineering, COMSATS University Islamabad, Islamabad Campus, Pakistan saud [email protected], T:+92 519 049 251

Department of Electrical and Computer Engineering, COMSATS University Islamabad, Islamabad Campus, Pakistan tariq [email protected], T:+92 519 049 184

Muhammad Arsalan

AN US

Division of Electronics and Electrical Engineering, Dongguk University, 30 Pildong-ro, 1-gil, Jung-gu, Seoul 100-715, Korea [email protected]

Muhammad A. Khan

School of Computing and Communications, Lancaster University, Lancaster, UK. [email protected]

M. A. Khalil

ED

M

Department of Computer Engineering,University of Lahore, Lahore, Pakistan [email protected]

Abstract

AC

CE

PT

Glaucoma along with diabetic retinopathy is a major cause of vision blindness and is projected to affect over 80 million people by 2020. Recently, expert systems have matched human performance in disease diagnosis and proven to be highly useful in assisting medical experts in the diagnosis and detection of diseases. Hence, automated optic disc detection through intelligent systems is vital for early diagnosis and detection of Glaucoma. This paper presents a multi-parametric optic disk detection and localization method for retinal fundus images using region-based statistical and textural features. Highly discriminative features are selected based on the mutual information criterion and a comparative analysis of four benchmark classifiers: Support Vector Machine, Random Forest (RF), AdaBoost and RusBoost is presented. The results of the proposed RF classifier based pipeline demonstrate its highly competitive performance (accuracies of 0.993, 0.988 and 0.993 on the DRIONS, MESSIDOR and ONHSD databases) with the state-of-the-art, thus making it a suitable candidate for patient management systems for early diagnosis of the Glaucoma. Keywords: AdaBoostM1, Glaucoma, RusBoost, Random forest, Support vector machine

1

1. Introduction

7 8

2 3 4 5 6

Nowadays, people are unaware of the visual impairment and 9 blindness lesions of an eye that are macular degeneration, hy- 10 pertension, glaucoma, and diabetic retinopathy (Soomro et al., 11 2018; Lee et al., 2015). Glaucoma is an incurable eye pathology 12 around the optic disc (OD) that badly damages the optic nerves 13 14

✩ An

Explicit Method for Localization of Optic disc from Retinal fundus

Images Preprint submitted to Expert Systems with Applications

15 16

and leads to vision loss. Once started, it cannot be cured but can be prevented by diagnosing at its early stage (Soomro et al., 2017b,a). According to an estimate, it becomes the most crucial and leading cause of vision blindness, it may affect about 80 million people by 2020 (Quigley & Broman, 2006; Khan et al., 2017, 2018). The progression rate of this disease is slow but occurs gradually over an extended time interval. The symptoms of glaucoma appear after an extensive time period when the disease is quite advanced. It is incurable but the growth rate of its progression rate have been regressed by proper treatment December 4, 2018

ACCEPTED MANUSCRIPT

23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61

Figure 2: Typical challenge in OD localization from retinal fundus. Green circular lines: experts mark OD boundary; blue circular lines: represents the OD by (Yin et al., 2011) and (Cheng et al., 2011) in both of the examples; skyblue lines mark disease around OD, where one of them is the PPA boundary.

To enhance the OD region and suppress the vessel background, appropriate image preprocessing is applied to increase the generalization ability of the classifier. The proposed method utilizes the statistical and textural properties of regions to learn a discriminative classifier from positive and negative regional samples. Pixel-accurate ground truth annotations are employed to generate training instances. In this work, various discriminative classifiers are trained and the Random Forest classifier is selected as the proposed approach due to its better generalization performance. The proposed method can robustly localize OD with high accuracy, which makes it a suitable candidate for glaucoma screening and diagnosis. The specific contributions of this work are:

AN US

22

M

21

ED

20

PT

19

if it is diagnosed at its primary stage. The glaucoma progression is symptomized from the few signs of vision loss. In Australia, 50% of peoples are unfamiliar with glaucoma (Maldhure & Dixit, 2015). The OD is defined by its central bright region called optic cup (OC) and the corresponding outer peripheral region called the neuroretinal rim, as shown in Fig. 1 . The cup to disc ratio (CDR), which is the ratio of vertical cup diameter 62 (VCD) to vertical disc diameter (VDD) is considered as a vi63 tal factor in glaucoma screening (Almazroa et al., 2015), due 64 to its directly proportional relationship with glaucoma occur65 rence. The segmentation of the OC has been demonstrated to 66 rely on accurate localization of the OD in literature (Yin et al., 67 2012). On similar lines, in this work, we argue that the correct 68 localization of the OD is vital for eye pathology screening and 69 diagnosis. 70 Expert and intelligent systems continue to contribute towards 71 scientific and technological advances in numerous fields from 72 medical applications, telecommunications, economics, trans- 73 portation, and surveillance. An important milestone for intel- 74 ligent machine vision systems is to match or surpass human vision performance. Expert automated vision systems have been successfully applied to disease diagnosis: recent works 75 include real-time detection of tuberculosis using an intelligent 76 mobile-enabled expert system (Shabut et al., 2018), deep neural 77 network based recommender system for skin lesion extraction 78 (Soudani & Barhoumi, 2019) and a convolutional neural net- 79 work based segmentation framework for breast tumor classifi- 80 cation (Rouhi et al., 2015). 81 The multifaceted challenges in OD localization include deformable shape, variation in color, OD boundary discontinuities 82 due to blood vessels and glaucoma pathology such as peripap- 83 illary atrophy (PPA) (Jonas et al., 1992), and ISNT rule (Hariz- 84 man et al., 2006). Most of the unsupervised and supervised 85 methods for OD boundary detection are gradient dependent ap- 86 proaches. Hence, they struggle to capture true OD boundary in cases where peripapillary atrophy (PPA) makes the OD bound- 87 ary discontinuous or the vascular structure originating from the 88 OD could misguide the segmentation. Also, the supervised ap- 89 proaches for OD localization (Fan et al., 2018) typically give 90 less attention to important stages of the learning pipeline, in- 91 cluding preprocessing, appropriate feature selection and data 92 imbalance problem in medical images. 93 To this end, we proposed a supervised approach for OD lo- 94 calization which is independent of boundary information, the 95 deformable shape of the OD and the irregularities in its size. 96

CE

18

AC

17

CR IP T

Figure 1: Optic disc and optic cup main structures. The enclosed region with green circle is for optic disc, central region enclosed with blue circular is the optic cup.

3

1. The problem of OD segmentation is modeled as a regional classification problem through a systematic classification pipeline. 2. Image regions are characterized by statistical as well as textural properties in a well-formulated classification framework to make it robust against multifaceted challenges in OD localization. 3. A thorough comparison of regional classifiers is presented for OD localization in retinal fundus images. Our results demonstrate the improved generalization of the Random Forest classifier over other classification based methods and unsupervised counterparts. The rest of this paper is organized as follows. A detailed review of the relevant methods is presented in section 2. The proposed method’s main steps for OD localization are presented in sections 3-8 that comprise of the preprocessing, superpixel segmentation, feature extraction and selection and the classification of each superpixel by using machine learning tools. Sections 9 and 10 comprise of experimental setup that contains datasets description, parameter selection, and comparative experimental results. The discussion and conclusions are described in sections 11 and 12, respectively.

CR IP T

ACCEPTED MANUSCRIPT

Figure 3: The train and test phases of the proposed method. The preprocessing and feature engineering blocks are common to both the train and test phases.

97

2. Literature review

138 139

105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137

A morphology-based method was proposed by (Walter et al., 2002) that extracts OD using watershed transformation. This method is based on a superfluous assumption that the OD is usually represented by bright pixels in fundus images, which stands invalid for pathology and real-life screening images. In (Morales et al., 2013), the gray image obtained by PCA is chosen as the input. The stochastic watershed is applied to extract the watershed regions. Afterward, region discrimination is performed to select the pixels which belong to OD based on the average intensity of the region. The process of region discrimination on the basis of the average intensity of regions in the watershed segmentation step could fail in scenarios where OD boundaries are not well-defined, as in case of Glaucoma affected images. Also, the initialization of the markers is crucial for the efficient working of the approach proposed by (Morales et al., 2013). (Welfer et al., 2010) proposed an adaptive method for the segmentation of the OD using an adaptive morphological approach. They used the particular attributes of the OD to address the issue including solid distracters along the vessel edges and peripapillary decay. Although, the approach was promising, the overall performance was disproportionate due to the deformable arbitrary shape OD regions obtained as the output of the watershed segmentation step.

AN US

104

M

103

ED

102

PT

101

CE

100

Computerized OD localization and segmentation is still con-140 sidered as a highly challenging task and an open problem in medical image analysis and diagnosis. The obvious challenges141 include the ever-varying physical characteristics of OD (such142 as size, structure, vessel structure, color) and various glaucoma143 pathology, as shown in Fig. 2. Considering only the indepen-144 dent scanning protocols, all of the mentioned properties are un-145 known. Many OD and eye lesion detection algorithms have146 147 been developed for 2D retinal fundus images. The current state-of-the-art techniques on OD localization148 and segmentation can be categorized into two categories: un-149 supervised and supervised approaches. (Morales et al., 2013)150 reviewed the unsupervised based OD segmentation methods in151 three ways, namely template based (Park et al., 2006; Lalonde152 et al., 2001), deformable model based (Morales et al., 2013;153 Lalonde et al., 2001; Lowell et al., 2004), and morphology-154 based (Morales et al., 2013; Walter et al., 2002; Welfer et al.,155 2010). All these methods finally mark the disc boundary by cir-156 cular Hough transforms due to its computational efficiency. The157 circular ellipse is fitted for glaucoma screening. However, all158 these methods are unsupervised. Some other methods that are159 pixel classification based methods (Welfer et al., 2010), (Mani-160 nis et al., 2016b) and hybrid soft computing based methods are161 162 also used for OD segmentation. 163 In the template-based matching methods, (Park et al., 2006) proposed the OD segmentation by thresholding followed by the164 circular Hough transform. The region with the highest average165 intensity indicates the OD. The threshold-based approach in-166 volves many parameters including the thresholds, which cannot167 be suited to all images. Also, the highest intensity assump-168 tion can be misleading in pathology images or in screening.169 (Lalonde et al., 2001) adopted canny edge detector to extract170 the boundaries of OD and Hausdorff template matching is used171 to fit the OD boundary by a circle. Although, the results were172 promising, however, the suitability of the thresholds for the173 edge detection process for a wide variety of images is question-174 able. (Lowell et al., 2004) proposed a deformable model-based175 method for OD segmentation by employing the image gradient176 based optimal points on the active contour of the image edges.177 Although, it was the first most comprehensive method for OD178

AC

98 99

localization and segmentation, however, the deformable model is sensitive to noise due to peripapillary atrophy and pathologies.

4

A supervised method for OD segmentation based on structured learning is proposed recently by (Fan et al., 2018). This method trains a structured forest for OD edge detection that is based on the Random Forest classifier. The OD edges obtained by the trained structured forest are thresholded, postprocessed and a circular Hough transform is applied to obtain the OD boundary. This is a competitive approach for OD detection and achieved performance comparable to the state-ofthe-art methods on benchmark datasets. However, the empirically sought parameters of the thresholding step and the circular Hough based ellipse fitting stage limit its generalization to unseen images. Also, it would not be the best of methods for cases in which the optic disc do not have a distinct edge due to disc tilting or peripapillary atrophy (PPA) and disc vessels could misguide the segmentation process .

ACCEPTED MANUSCRIPT

181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200

Recently deep neural network approaches have been inves-221 tigated for optic disc segmentation (Maninis et al., 2016a; Fu et al., 2018; dos Santos Ferreira et al., 2018) with exceptional222 performance on both separate-training and cross-training eval-223 uations. Deep neural networks based approaches are highly224 promising for real-time screening and diagnosis due to their225 high accuracy and generalization on unseen data. However, this high accuracy is achieved at the cost of huge training data and226 high computational overhead of training. The generation of a227 large set of training examples involves systematic data augmen-228 tation. However, the challenge lies in generating variations that229 capture the underlying distribution of a high degree of varying230 fundus images. The high degree of inter dataset variation makes231 this problem highly challenging in scenarios where there is only232 a limited number of available training examples. The com-233 putational requirements of the training process of deep neural234 networks is yet another factor which limits their application in235 screening and diagnosis applications. While large clinical facil-236 ities could house high-end computing resources; not all screen-237 ing facilities have such powerful computing capabilities. Also,238 the large training times make simultaneous training from new239 patients data impractical. 240 241

206 207 208 209 210 211 212

213 214 215

216 217 218 219 220

3.1. Preprocessing In image processing, preprocessing is a preliminary step that is used to enhance the image according to the requirement for feature extraction. In this paper, this has been subdivided into three subsections. Our retinal fundus image preprocessing proposed model includes noise removal, image enhancement, and image cropping. The first step of preprocessing is to smooth the image by removing the noise and enhancing the edges of the optic disc. For this purpose, bilateral filtering is used. A bilateral filter is a non-linear, edge-preserving and noise smoothing filter. It replaces the intensity of every pixel with a weighted normal of intensity values from adjacent pixels. This weight is founded on a Gaussian approximation. The bilateral filter is likewise characterized as a weighted normal of adjacent pixels, in a way fundamentally the same as Gaussian convolution. The difference is that the bilateral filter considers the distinction in values with the neighbors to protect edges while smoothing. The formalization of the bilateral filter is given in Eq. 1

M

205

ED

204

PT

203

For offline diagnosis and annotation, where there is no time constraint, deep neural networks or ensemble methods would251 result in higher accuracy as compared to the proposed method. 252

CE

202

The proposed method uses DRIONS and also further evaluated on MESSIDOR and ONHSD dataset. The main block diagram of the proposed work is depicted in Fig. 3, the details of these block are described in the following subsections.

The proposed method is favorable in scenarios where there242 is scarcity of training data. As the proposed method is region-243 based, limited number of training images could suffice to gener-244 ate a reasonable number of training instances. The low training 1 X time and minimal computational requirements make the pro-245 Gσs (kp − nk)Gσr ( Ip − In )In , (1) Np (I) = wp n∈W posed method a suitable candidate for online training and updating on new patients data. Another noteworthy strength of the proposed method is that unlike previous approaches, it uti-246 where  Gσs (kp  − nk) is a gray level similarity function, lizes multiple cues and modalities to obtain a region-wise pre-247 Gσr Ip − In is a geometric closeness function, Np is the fildiction. This essentially enables it to distinguish OD from other248 tered image and wp controls the normalization of the filtered unhealthy retinal structures and provides robustness to pathol-249 output, given as ogy images. X   wp = Gσs (kp − nk) Gσr Ip − In . 250 (2) 253

We note that the accuracy of the proposed regional classifi-254 cation pipeline is dependent upon the OD localization ability of255 the proposed approach. If the localization scheme fail on some256 images, the regional classification pipeline is highly likely to257 258 obtain sub-optimal results.

AC

201

3. Method

CR IP T

180

AN US

179

n∈W

The vectors p and n are the pixel locations being considered in the window W, where Ip and In are the image intensities at these pixel locations, respectively. The next step in preprocessing is image localization and cropping around the OD region. Localization is accomplished by finding the centroid of the most intense and least eccentric region of the green channel preprocessed image. The localization process of the OD is generally invariant to the color channel

Figure 4: Preprocessing. (a) Original image cropping (b) Original image with Red channel (c) Good contrast image for reference (d) Reference image with Red channel (e) Final preprocessed histogram matched image.

5

ACCEPTED MANUSCRIPT

276

315

266 267 268 269 270 271 272 273 274

277

278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298

Ccrop =

M N 2 XX I(x, y), M × N x=1 y=1

(3)

316

3.2. Superpixel Segmentation The image is segmented into superpixel by using Simple Linear Iterative Clustering (SLIC) technique. This technique is proposed by (Achanta et al., 2012) for partitioning the image into equally sized square shape. In SLIC, some parameters are required which can be selected by visualizing the variations between them and boundary observance. For initial superpixels, R is considered as superpixel region size. The center coordinates of each superpixel are iteratively updated. The grouping of pixels are based on their distance metrics of spatial and intensity values of distances. Spatial distance D s between two consecutive pixel (e.g ath and bth )is calculated in Eq. 4 q (4) D s = (xa − xb )2 + (ya − yb )2 ,

where a and b are the pixel locality components. The intensitybased distance Di between the two adjacent pixels is defined in Eq. 5 p Di = (Ia − Ib )2 , (5)

AN US

265

317

318

where x and y are the pixel coordinates and M and N are the319 width and height of the image. The thresholding criterion of Eq. 3 is widely adopted in the object detection literature320 (Achanta et al., 2009) to segment important objects and sup-321 press background information. This adaptive threshold is ex-322 pected to only retain the high-intensity objects such as the optic323 disc or pathologies and suppress all the unwanted background324 noise. This adaptive threshold is found to be robust in our experiments as it does not have any negative effect on the eccen-325 tricity computation step. Centroid of all regions, obtained after thresholding, are calculated. Eccentricities of these regions are326 compared and the centroid of the one with the lowest eccentric327 value is used for cropping. Subimage size is experimentally328 set to be 200×200. Multiple windows sizes are placed over the329 localized centroid and 200×200 stands optimum because it nei-330 ther catches extra information nor it trims the OD edges. 331 The last preprocessing step is color channel selection and332 histogram matching. In color image, empirically, three stan-333 dard channel of color are considered mostly, i.e. (Red, Green,334 Blue). A color image is formed by the combination of these335 three channels. Most of the information for OD is in the Red336

M

264

ED

263

PT

262

CE

261

channel image, therefore, this channel is used for further processing. Background information image is normalized by subtraction with its estimated background. After that we have applied the histogram matching (Gonzalez, 2006) to normalize the image variations.

where Ia and Ib are the normalized intensity values of the ath and the bth pixel, respectively. The total distance measure, which is the combination of spatial and intensity distances, is then calculated in Eq. 6 r Ds (6) D = Di + ( )2 r2 , R where r is a regularizer coefficient. A regularizer coefficient defines the flexibility of superpixel boundaries. If the value of r is higher, it results in more regularize segments. A lower value of r creates boundaries that are more flexible. For obtaining optimal regularizer coefficient r, the intensity values of retinal images in Eq. 6 are normalized in range [0, 1]. The normalization is performed to ensure that both of the distances (i.e. spatial and intensity) are in the same range. Fig. 5 shows cropped retinal image containing an optic disc, which is segmented into small patches of two different region sizes, R = 5 and R = 10. For both sizes, we used regularizer factor r = 0.2.

AC

260

CR IP T

275

of the preprocessed image. But in some cases, the performance299 of green channel for localizing OD is better than Red channel.300 Red channel usually posses the information of contrast, in other301 words, high contrast regions are more prominent in Red region.302 In case of retinal vessels image, the high contrast regions (optic303 disc and pathologies) are more prominent in Red channel. In images where pathologies are dominating, this channels fails to304 localize the true OD region as it possess pathologies having less305 eccentricity than OD. In contrast, in green channel, the OD of306 region is more eccentric than pathologies, for such cases. Ec-307 centricity is the measure of a circular shape, higher the region308 shape closer to circle, lower the value of eccentricity. For this309 purpose, the green channel is morphologically dilated with a310 disk because OD is almost disk-shaped. Disk size is empiri-311 cally set to be 20. A threshold computed as twice the mean312 intensity of the image is applied to the image to get a binary313 image. The thresholding criterion is adaptive which is given by314

259

Figure 5: Example over-segmentations with various window sizes: (a) retinal fundus image, (b) superpixel with region size R=5 (initial grids of 5 × 5), (c) superpixel with R=8, (initial grids of 8 × 8) with fixed r=0.2

6

ACCEPTED MANUSCRIPT

Figure 6: Intensity based statistical features. (a) Average intensity feature, (b) Maximum intensity feature, (c) Minimum intensity feature, (d) Median intensity feature, (e) Mode intensity feature

4. Feature Extraction

361 362

340

Feature extraction and normalization steps are very important363 to make the classification based approach more robust for OD 364 segmentation that are detailed below.

341

4.1. Intensity based statistical features

338 339

365

349

4.2. Texton-map Histogram

346 347

350 351 352 353 354 355 356 357

372 373

Retinal images have a very complex OD structure. For this374 reason, we do not rely solely on intensity characteristics and therefore additionally employ texton-map based histogram fea-375 tures to make OD segmentation more robust. The texton-map376 basically capture the texture of images. The textural informa-377 tion is extracted from the image by convolving it with a bank378 of filters. The filter bank used here is created from Gabor filter.379 380 Gabor formulation is (Henriksen, 2007) defined in Eq. 7:

M

345

ED

344

PT

343

0

360

0

381

384

The σ defines the size of filter, wavelength of the sinusoidal 385 factor is denoted by λ, phase shift is ψ, and γ for aspect ratio

CE

359

0

−a 2 + γ2 b 2 a 382 G(a, b; θ, λ, ψ, γ) = exp( ) exp(i(2×π +ψ)).(7) 2 383 λ 2σ

AC

358

0

a = a cos θ + b sin θ, 0

b = a sin θ + b cos θ.

(8) (9)

The parameter used to tune the Gabor filter bank (GFB) are discussed in more detail in the section of texton-map parameters. The retinal fundus images are convolved with the kernels of GFB. The procedure for GFB is given in Fig. 7. After convolution, response vector for all filters in the bank is generated. This response vector is then clustered into k clusters by using k-means. Features are extracted from each superpixel by the histogram value of texton-map for that superpixel.

AN US

348

Intensity based statistical features belong to the class of first 366 order statistical features (Jain, 1989). These features directly 367 deal with the gray level of images within the specific region 368 (i.e. region of interest (ROI)), that are superpixels in this work. 369 Five intensity features: average, maximum, minimum, median 370 and mode intensity, are extracted for each superpixel. The vi371 sualization of these features is presented in figure 6.

342

of spatial components. In Eq. 7, θ represents filter kernel direc0 0 tions that are used to calculate the terms a and b defined in Eq. 8 and in Eq. 9:

CR IP T

337

386 387 388 389 390 391

392

393 394 395 396 397

Figure 7: Process of the texton map computation based on the Gabor filter bank.398

7

4.3. Fractal features

The fractal features are extracted by the binary decomposition of images, using a multilevel Otsu thresholding algorithm. It is important to select the optimum value of the threshold. The desired value of thresholds is discussed in parameter selection section. Otsu (Liao et al., 2001) is a multilevel thresholding method and these levels depend on the threshold value. After thresholding the image, edge detection is applied to extract the boundaries. The spatial representation for fractal features is depicted in Fig. 8. From each binary Otsu image and fractal border image, three characteristics are extracted: average intensity, area and fractal dimension, (Canny, 1987). The area is the sum of edge pixels in each superpixel, the mean intensity is the mean value of each superpixel in the entire threshold region. The fractal distances or Hou’s dimensions are used to find the complex binary region within each superpixel, calculated by plotting the bounding box on boundaries. The box formulation is given in Eq. 10. =0 = lim

→0

log n() , log  −1

(10)

where n() denotes the number of the box with side length represented by . An approximation of fractal dimension is done by using a box counting algorithm. Fig. 9 shows the example of fractal features, extracted from the binary Otsu decomposition and canny edge image for each superpixel of the retinal fundus image.

CR IP T

ACCEPTED MANUSCRIPT

AN US

Figure 8: The flow of fractal feature computation.

403 404 405 406 407 408

409

410

411 412 413 414 415 416 417 418 419 420 421

ED

402

PT

401

For classification-based approaches, feature normalization (Cao et al., 2016) is important because most of the classifier work on distance based schemes. It is necessary to obtain a robust classification by normalizing the feature variable in a specific range. There are many approaches that are used to standardize the features by mean and histogram normalization. In this paper, mean normalization is used to standardize the features. Given feature f , the mathematical formulation for mean 422 normalization is given in Eq. 11 and 12. Meannorm =

CE

400

5. Feature Normalization

f − fmin , fmax − fmin

AC

399

M

Figure 9: Otsu and Canny based fractal features.

Normalization = Meannorm × (R − D) + R.

(11) (12)

423 424 425 426

In Eq. 12 R denotes the maximum of the range of the desired427 target, whereas D denotes the minimum of the range of the de-428 sired target. fmax and fmin represent the maximum and mini-429 mum values of the feature f . In summary, total 17 features are extracted from each super-430 pixel. Feature matrix includes 5 intensity statistical, 6 textonmap features and 6 fractal feature (i.e. 3 from each binary431 decomposed image). It is important here that the range value432 for normalization [0, range(texton-map)] is selected from the433 texton-map feature. In this paper, it is empirically selected by434 435 visualization. These are also depicted in table 1. 8

Table 1: Total number of extracted features

Feature Name Statistical based features Texton histogram Fractal features Total

Number of Features 5 6 6 17

6. Feature selection Feature selection is an essential step for supervised methods. Feature selection helps in fighting with the curse of dimensionality. Good features improve the prediction performance of the classifier. In this paper, features are selected on the basis of features mutual information (MI(a,b) ) to find the minimum redundancy between feature sets (Peng et al., 2005). Formally, the MI is defined as follows: XX ℵ(a, b) , (13) MI (A, B) = ℵ(a, b). log ℵ(a.b) a∈A b∈B where ℵ(a, b) is the joint mass probability of feature a and b and ℵ(a, b) = ℵ(a.b) when a and b are statistically independent. Features are discarded based on the minimum redundancy criterion of (Peng et al., 2005) until the minimum number of selected features is reached.

ACCEPTED MANUSCRIPT

7. Classification

477

478 479

443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476

The image segmented by the classifier is binarized for further processing. As vessels are also present in the OD region, they also affect the OD region during binarization. In addition, the boundaries of OD can be affected by the nonuniform illumination. To mitigate these effects and estimate the true OD region, the largest connected object is obtained and its boundary is used as the raw estimation. The best-fitted ellipse is computed as the OD boundary (Zhang et al., 2009). The ellipse fitting minimized the noise effects introduced by vessels especially from the inferior and superior regions of the OD neuroretinal rim. Fig. 10 shows the effects of OD location by raw estimated boundary in Fig. 10 (b) and ellipse fitted boundary in Fig. 10 (c). 9. Experimental analysis

In this section, the results of four nonlinear classifiers are evaluated on the publicly available retinal fundus datasets DRIONS, ONHSD and MESSIDOR to test the robustness of the proposed method. This section reports our experiments including parameter selection and results evaluation.

AN US

442

M

441

ED

440

PT

439

In supervised classification the model is constructed by learn-480 ing from data along with its annotated labels. The learned481 model is then evaluated on unseen data. In the proposed work,482 we employ four supervised classification methods for optic483 disc segmentation and present a comprehensive comparison on484 benchmark datasets. 485 The Support Vector Machine (SVM) classifier is a supervised486 classifier used for classification (Boser et al., 1992). It can be487 used for both linear and nonlinear classification by using its488 kernels. SVM generates a hyperplane in a high-dimensional489 feature space, where the objective is to maximize the hyper-490 plane margin from the support vectors. In this work, a radial basis function is adopted as the kernel, due to the nature of the 491 problem. The random forest (RF) classifier constructs multiple deci-492 sion trees based on the input variables. The RF classifier op-493 erates according to the principle of bootstrap aggregation (i.e.494 bagging) to increase the accuracy of generalization and reduce495 overfitting. Bootstrap sampling is applied to the construction496 of trees where a number of random predictors are utilized at each decision split. The idea is to add the responses of multiple497 trained trees in such a way that the correlation between trees 498 is increased, while the variation between trees is reduced. At 499 the test stage, the responses from multiple trees are aggregated 500 using majority voting (Breiman, 2001). RF has become an im501 portant data analysis tool in a short period of time, and became 502 popular because it can be applied to non-linear and higher order 503 data sets (Strobl et al., 2007). 504 The AdaBoost (Freund & Schapire, 1997) classifier aggre-505 gates the response of several weak learners to increase the gen-506 eralization accuracy. The notable difference is that they are507 reweighted instead of sampling (Breiman, 2001). In addition, a508 weak learner could be compared to an RF decision tree. 509 Random undersampling (RusBoost) is suitable for classify-510 ing imbalanced data when the instances of one class dominate511 many times more than the other. Machine learning techniques512 fail to efficiently classify skewed data, but RusBoost solved the513 problem by combining sampling and boosting. We explored514 these classification algorithms, and the results are reported in515 the result section. 516

CE

438

AC

437

8. Ellipse Fitting

CR IP T

436

9.1. Dataset description The DRIONS dataset contains 110 color retinal fundus images and annotation of optical nerve head (ONH) by two trained human experts. The data was obtained from an analogical color fundus camera and digitized using a high-resolution scanner HP-Photo Smart-S20. The size of each image is 8 bits per pixel with a resolution of 600×400. The average age of the patients in this database is 53.0 years with a male to female ratio is 46.2%. The ratio of patients with chronic glaucoma disease is 23.1%. The manual annotated OD segmented by the first expert has been considered as a gold standard in this paper. The MESSIDOR data set contains 1200 retinal images that were collected from three different ophthalmological departments. Out of 1200 images, 800 were obtained with a dilatation of the pupil with tropicamide at 0.5 % and were captured with the 3CCD Topcon TRC camera. ONHSD database contains 99 retinal fundus images that are taken randomly from 50 patients with diabetic retinopathy; out of 99 images 96 are discernable ONH. The resolution of images is 640×480.

Figure 10: The procedure of ellipse fitting. (a) Original retinal image, (b) Classifier prediction overlayed on original image, (c) Ellipse fitting.

9

CR IP T

ACCEPTED MANUSCRIPT

M

AN US

Figure 11: Representatives retinal fundus images with different types of light artifacts in upper row and their ground truth manual segmentation in the lower row.

ED

Figure 12: Superpixel based over-segmentation with R = 8, where several variations based on different values of the regularization parameter are: (a) r = 0, (b) r = 0.2 and (c) r = 0.5.

517 518

PT

519

Tree Depth Evaluation 0.15

521 522

CE

0.13 0.12 0.11

523 524 525 526

AC

Out−of−Bag Classification Error

0.14

520

527

0.1

528 529

0.09

530

9.2. Parameter selection The determination of parameters is an essential step in the classification-based learning problems. In feature extraction, some features are nonparametric such as statistical features that are computed directly from intensity values of superpixels. The parameter for superpixel segmentation, texton map, and fractal features are imperative. In SLIC, the size and regularization parameters control the shape and size of the superpixels. They must be appropriately set such that superpixels represent the object boundaries. In the calculation of texton feature, the parameters of the Gabor filter bank and the number of K-means clusters must be determined in an appropriate manner. In our proposed method, these parameters are set during the training phase.

0.08

531

0.07

532

0.06 0

20

40 60 Number of Grown Trees

80

100

533 534 535

Figure 13: Optimal tree depth selection based on minimal OoBError. The graph shows the relationship between the OoBError and the number of trees, where536 the optimal value for number of trees is observed to be 25, after which the error537 goes into steady state. 538

10

9.2.1. Superpixel parameters Important parameters for superpixel segmentation are the region size R and the regularizer r. The regularizer r = 0.2 is determined to be the best-suited value in our experiments, which results in superpixel shapes that represent OD boundaries, as shown in Fig. 12 . The region size R is determined empirically as a function of dice score and computational complexity in this work. Table 2 presents the dice score for various region sizes.

ACCEPTED MANUSCRIPT

Table 2: Optimal superpixel size selection for maximizing the Dice coefficient.

Size of Superpixel Dice coefficient

539 540

4

6

8

10

15

0.94

0.92

0.89

0.88

0.81

Table 3: Confusion Matrix.

OD Available Classify correctly Classify wrongly

R = 8 is chosen in our implementation as a trade-off between accuracy in terms of Dice score and computational complexity.

True Positive (TP) False Negative (FN)

OD Not Available False Positive (FP) True Negative (TN)

541

544 545 546 547 548 549 550 551 552 553 554 555 556

9.2.2. Texton Histograms parameters For the kernel direction in Gabor filter bank, six different ori-567 entations are selected: [0◦ , 30◦ , 45◦ , 60, 90◦ , 120◦ ]. These six568 directions are enough to cover the entire space with an appro-569 priate step size. By adding more directions we can get more570 information, but it also adds redundant information and effect571 computational complexity. The wavelength coefficient and fil-572 ter size are selected empirically in this work. The image with filter size 0.3 is close to the original image and appears to be573 quite blurred for values above 1.5. Therefore, filter sizes are574 selected in the range [0.3,1.5] with an increment of 0.3. The575 values of the wavelength coefficients that are selected by visual576 inspection of the filters are in the range of [0.8,1.5]. The num-577 ber of clusters ’k’ is set equal to the number of filter kernels,578 which is set to six in our implementation. 579 580

564

9.3. Out of bag error (OoBError)

562

ED

561

PT

560

CE

559

AC

558

M

563

9.2.3. Fractal feature parameters In the fractal feature computation, it is important to select581 an appropriate threshold level. In this paper nt = 3 is chosen 582 which creates 3 sub-levels of threshold. We selected this value owing to the considerably better Dice score for nt = 3. By583 increasing the levels of thresholds, the classification becomes584 more complex and time-consuming at fractal feature extraction.585

557

The prediction error of the classifier on the observations that are not used in the constructing the next base learner is known as the OoBError. OoBError is evaluated during the training stages of the learning classifiers, including RF, AdBoost and RusBoost. This empirical evaluation is sufficient to estimate the unbiased error and thus cross-validation is deemed unnecessary. Fig. 13 shows that the OoBError is minimum at 25 trees and the variation in OoBError is negligible afterwards. 9.4. Evaluation measures Five quantitative measures are employed in this work including sensitivity, specificity, Dice similarity coefficient, accuracy and area overlap. The quantitative measures are computed on the basis of the confusion matrix presented in Table 3. The quantities from the confusion matrix are explained below:

AN US

542 543

566

CR IP T

565

586

587

588 589 590

591

592 593 594

595

596 597

598

599

600

Figure 14: Visualization of the evaluation quantities in terms of predicted and ground truth regions. Table 3 gives better understanding of the evaluation mea-601 sures. 602

11

• TP: correctly identified OD pixels.

• FP: incorrectly identified non-OD pixels as OD pixels. • TN: correctly identified Non-OD pixels.

• FN: incorrectly identified OD pixels as non-OD pixels.

Dice coefficient is a similarity measure overlap ratio between ground truth (Bernal et al., 2017) and predicted output of the classifier (Zou et al., 2004). The mathematical formulation of dice coefficient is given in Eq.14. 2 ∗ TP (14) Dice = (2 ∗ T P + FN + FP) The true positive ratio is called sensitivity (SN) and also known as recall. It measures the portion of correctly identified TP as true class, given as TP SN = . (15) T P + FN Specificity (SP) is the true negative rate which means the number of negatives (healthy subjects) that are predicted as negatives. It can be measured as: TN SP = (16) T N + FP The classification accuracy is a sample of correctly identified observations in each class. TP + TN Accuracy = (17) T P + FP + T N + FN The area overlap (AOL) is defined as follows TP AOL = . (18) T P + FN + FP Ten-fold cross validation is employed for all classifiers in all of our experiments.

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 15: Visual results of the proposed method on representative images from the DRIONS database.

603

10. Hierarchical classification

623 624

610

10.1. Quantitative evaluation

611 612 613 614 615 616 617 618 619

620 621 622

ED

608

631 632

Table 5 gives a brief comparison of the regional classifiers633 on benchmark datasets. It can be observed that all the clas-634 sifiers obtain competitive results in terms of selected perfor-635 mance measures. In general, the average performance of the AdaBoostM1 classifier is better than that of the SVM in terms636 of AOL, Dice and Accuracy. Whereas the Random forest and637 the RusBoost classifiers obtain comparable performance no-638 tably higher than AdaBoostM1 and SVM. The major findings639 640 are summarized here:

PT

607

CE

606

AC

605

M

609

In this section we present the quantitative comparison of the625 classifiers for the task of OD segmentation in terms of the se-626 lected measures on all datasets followed by visual results of627 the proposed method on representative images from benchmark628 databases. Next, a comparison of the proposed method with the629 630 state-of-the-art methods is presented.

604

641

• The performance improvement of RusBoost over Ad-642 aBoostM1 can be attributed to its class imbalance handling643 644 before boosting.

• The Random forest classifier outperforms AdaBoostM1 on average in terms of all performance measures. This result is against the expectation as boosting generally performs better in practice than bagging. This result can be attributed to the fact that AdaBoostM1 learns from past misclassified examples, however, as OD segmentation is a class imbalanced problem, the AdaBoostM1 classifier slightly overfits to the majority class. • The average performance of the RusBoost classifier is comparable with the performance of the Random forest classifier. However, the low computational overhead of the Random forest classifier at the training stage makes it more feasible.

Table 4 presents a comparison of classifiers in terms of training time. Also, the individual times of the various stages of the proposed pipeline at the test stage are also reported in Table 4. It can be observed that the computation time at the test stage is dominated by the feature computation stage. The texton feature computation contributes the most to the overall feature computation time. Visual results of the proposed approach are presented in 15 and 16. The predicted boundaries by the classifiers are show

Table 4: Timing comparison of classifiers training time and stage-wise timing of the proposed classification pipeline.

Time (seconds)

Classifier Training

Feature Computation

Feature Selection + Classification

Support Vector Machine Random Forest AdaBoostM1 RusBoost

36.3 9.1 7.5 12.9

19.5

2.5

12

ACCEPTED MANUSCRIPT

646

AdaBoostM1

RusBoost

DRIONS

AOL Dice Acc SN SP

0.809 0.892 0.992 0.969 0.993

0.821 0.899 0.993 0.948 0.994

0.786 0.872 0.992 0.845 0.997

0.835 0.902 0.993 0.959 0.993

MESSIDOR

AOL Dice Acc SN SP

0.733 0.841 0.988 0.934 0.987

0.747 0.851 0.988 0.953 0.988

0.742 0.847 0.988 0.888 0.953

0.743 0.847 0.988 0.912 0.988

AOL Dice Acc SN SP

0.783 0.875 0.99 0.926 0.991

0.824 0.897 0.993 0.924 0.995

0.797 0.879 0.992 0.877 0.996

0.826 0.902 0.993 0.952 0.994

with different colors, while the ground truth boundary is shown652 653 in blue. 654

647 648

10.2. Comparison of performance measures with other algo-655 656 rithms 657

M

mark method selected for comparison include: the technique of (Fan et al., 2018), (Zahoor & Fraz, 2017),(Abdullah et al., 2016), (Morales et al., 2013) and (Walter et al., 2002). These method are chosen due to their state-of-the-art performance and recency. It is worth mentioning here that due to the highly competitive performance of the RF classifier and its relatively lower training time, all the following comparisons are based on choosing the Random forest as the regional classifier in our proposed

Table 6: Comparison of the proposed approach with the state-of-the-art methods in terms of accuracy.

1

Author

Method+Brief Detail

(Zahoor & Fraz, 2017)

2

(Fan et al., 2018)

3

(Walter 2002)

4

5

ED

No.

Morphology based preprocessing followed by circular Hough transform for OD localization

PT

651

CE

650

The performance of the proposed method is compared with 658 the benchmark methods in this section in terms of the five eval659 uation measures on three benchmark databases. The bench-

et

AC

649

CR IP T

Random forest

AN US

645

Support vector machine

ONHSD

Table 5: Comparison of RF, SVM, AdaBoostM1 and RusBoost classifiers for optic disc classification.

(Morales 2013)

et

Proposed

DRIONS MESSIDOR DRIVE

0.9986 0.9918 0.998

DRIONS MESSIDOR ONHSD

0.976 0.977 0.9895

OD extraction by watershed segmentation

DRIONS MESSIDOR

0.9689 —–

ONHSD al.,

Accuracy

Structured learning based OD segmentation

al.,

(Abdullah et al., 2016)

6

Database

PCA and mathematical morphology based OD segmentation

—–

DRIONS MESSIDOR

0.9934 0.9949

ONHSD

0.9941

Morphology based OD localization segmentation based on grow-cut algorithm

DRIONS MESSIDOR ONHSD

0.9909 0.9925 1

Multi-parametric optic disc Segmentation using superpixel based feature classification

DRIONS MESSIDOR ONHSD

0.993 0.998 0.997

13

ACCEPTED MANUSCRIPT

661 662 663 664 665

pipeline.

689

Table 6 presents the comparison of the proposed method with the state-of-the-art method in terms of classification accuracy.691 The results show that the proposed method obtains highly com-692 petitive results as compared with the state-of-the-art in terms of693 694 average accuracy for all datasets. 690

695

667 668 669 670 671 672 673 674 675 676 677 678 679 680 681

Table 7 presents a comprehensive comparison of the proposed approach with the state-of-the-art methods in terms of all selected performance measures. From Table 7 it can be ob-696 served that the sensitivity of the proposed method is consistently higher than all existing methods on all datasets. Also697 the accuracy of the proposed method is better than the state-of-698 the-art methods on average for all databases. The sensitivity of699 the proposed method is quite competitive with the state-of-the-700 art on all databases with measures of 0.994, 0.988 and 0.995701 on DRIONS, MESSIDOR and the ONHSD databases, respec-702 tively. The dice score of the proposed method for the DRIONS703 database is better than the methods of (Zahoor & Fraz, 2017)704 and (Walter et al., 2002), while comparable to the rest of the705 benchmark methods. Similarly, the dice score of the proposed706 method for the ONHSD dataset is comparable to the state-of-707 the-art methods and slightly better than (Morales et al., 2013).708

682

688

M

687

This paper presents a supervised approach for OD localization. The presented approach used DRIONS database, ONHSD and MESSIDOR data for evaluation and relies on intensity, texton-map histogram and fractal features. In the progression of feature computation process, we considered intensity, texton-map, and fractal features. Important parameters required at the feature computation stage were sought empirically as optimization of these parameters is beyond the scope of this work. Future work will focus on incorporating more discriminative features in the proposed method for regional optic cup classification. Accurate prediction of optic cup along with optic disc is necessary to obtain reliable cup-to-disc ratio, which is an important parameter in glaucoma screening. In superpixel based OD segmentation, the superpixel size and regularization are important parameters that directly effect the accuracy and computational complexity of the proposed classification pipeline. Furthermore, the scale and shape of the superpixel can lead to inaccurate label assignment on a regional level

Table 7: Comparison of the proposed method with the state-of-the-art methods in terms of Sensitivity, Specificity, AOL, Dice and Accuracy measures on three benchmark datasets.

ED

686

Performance measures→ Method↓

PT

685

Sensitivity

Specificity

AOL

Dice

0.938 0.8957 0.6715 0.928 0.8508 0.969

Accuracy

0.999 —– —– —– 0.9966 0.994

0.884 0.8473 0.6227 0.842 0.851 0.821

0.889 0.9137 0.6813 0.908 0.9102 0.899

0.998 0.976 0.9689 0.993 0.9549 0.993

0.889 0.9212 0.93 0.8954 0.948

0.997 —– —– 0.9995 0.988

0.844 0.8636 0.8228 0.879 0.747

0.903 0.9196 0.895 0.9339 0.851

0.991 0.977 0.9949 0.9989 0.988

0.9077 0.931 0.8857 0.924

—– —– 0.9992 0.995

0.8346 0.8045 0.861 0.824

0.9032 0.8867 0.9197 0.897

0.9895 0.9941 0.9967 0.993

DRIONS

(Zahoor & Fraz, 2017) (Fan et al., 2018) (Walter et al., 2002) (Morales et al., 2013) (Abdullah et al., 2016) Proposed

CE

684

709

The overall results of the proposed method are quite promis-710 ing as compared with the state-of-the-art. This improved per-711 formance of the proposed approach can be attributed to the fact712 that the proposed method is independent of a particular feature713 modality and simultaneously considers multiple OD attributes714 in a classification framework on a regional basis, thus improv-715

AC

683

11. Discussion

AN US

666

ing its generalization on unseen fundus images. In contrast, majority of the morphology and the segmentation based approaches depend upon the assumption of OD being the brightest region, thus do not generalize well on pathology images. Similarly, the supervised method of (Fan et al., 2018). is solely dependent upon edge information of the OD, which is an impractical assumption for peripapillary atrophy images.

CR IP T

660

Messidor (Zahoor & Fraz, 2017) (Fan et al., 2018) (Morales et al., 2013) (Abdullah et al., 2016) Proposed ONHSD (Zahoor & Fraz, 2017) (Morales et al., 2013) (Abdullah et al., 2016) Proposed

14

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 16: Visual results of the proposed approach for images from the MESSIDOR and ONHSD databases. The first row presents the visual results for SVM, RF, AdaBoostM1 and RusBoost classifiers on the MESSIDOR dataset, whereas the second row presents the results for ONHSD database.

721 722 723 724 725 726 727 728 729 730 731

732

733 734 735 736 737 738 739 740 741 742 743 744 745

M

720

ED

719

PT

718

and can cause the superpixels to miss OD boundary. Searching746 for an optimal region size and regularization parameters is an747 interesting problem, however, beyond the scope of this work. In748 this work, suitable values for these parameters are determined749 empirically. 750 We utilized only the five extracted features that have been751 chosen on the basis of their mutual information with four dif-752 ferent classifiers (SVM, AdaBoostM1, RusBoost, and RF) as in753 Table 5 and observed that the Random forest classifier offers the754 best trade-off in terms of high accuracy and low computational755 complexity at the training stage. A comprehensive comparison756 of the proposed framework with the benchmark methods (in757 terms of standard performance measures on widely accepted758 databases) demonstrate that its performance is comparable to759 the best performing benchmark methods and superior to several state-of-the-art methods. 12. Conclusion

CE

717

760

761

AC

716

This paper presented a regional classification framework for762 accurate localization and segmentation of the optic disc. The763 modeling of the optic disc segmentation as a region-based clas-764 sification by utilizing multi-modality attributes proved to be ro-765 bust against the multifaceted challenges of optic disc detection.766 The results demonstrate that the proposed method is resilient767 against the highly varying nature of optic disc appearance as768 compared with other expert and intelligent systems, which fail due to their reliance on a single feature modality. The proposed method obtained average improvement of769 around 0.3% upon the performance of the best performing state-770 of-the-art method and around 1.7% against the second best per-771 forming benchmark method in terms of accuracy. The notewor-772 15

thy improvements in terms of sensitivity and accuracy obtained by the proposed method as compared with the state-of-the-art expert methods can be attributed to its independence to a particular feature of the OD. To incorporate simultaneous optic cup detection for Glaucoma screening, the proposed classification framework is to be extended to multiclass classification framework in future by incorporating useful regional features for effective discrimination of cup regions and a mutilabel classification loss formulation. Another avenue to explore in future is the class imbalance handling for both optic disc and the optic cup regions. As the feature computation stages are independent, the parallel implementation of the proposed method for obtaining real-time performance is also a subject of our future research. 13. Acknowledgment The authors would like to specially thank the management teams of the MESSIDOR, DRIONS and ONHSD databases for providing retinal fundus image databases along with the manual expert annotations. These database were download from there publicly available links. The authors would like to thank Higher Education Commission Pakistan for funding the research reported in this paper under grant number 212020/SRGP/R&D/HEC/2018. References Abdullah, M., Fraz, M. M., & Barman, S. A. (2016). Localization and segmentation of optic disc in retinal images using circular hough transform and grow-cut algorithm. PeerJ, 4, e2003.

ACCEPTED MANUSCRIPT

781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843

CR IP T

780

retinal image understanding. In International Conference on Medical Image Computing and Computer-Assisted Intervention (pp. 140–148). Springer. Morales, S., Naranjo, V., Angulo, J., & Alca˜niz, M. (2013). Automatic detection of optic disc based on pca and mathematical morphology. IEEE transactions on medical imaging, 32, 786–796. Park, M., Jin, J. S., & Luo, S. (2006). Locating the optic disc in retinal images. In Computer Graphics, Imaging and Visualisation, 2006 International Conference on (pp. 141–145). IEEE. Peng, H., Long, F., & Ding, C. (2005). Feature selection based on mutual information criteria of max-dependency, max-relevance, and min-redundancy. IEEE Transactions on pattern analysis and machine intelligence, 27, 1226– 1238. Quigley, H. A., & Broman, A. T. (2006). The number of people with glaucoma worldwide in 2010 and 2020. British journal of ophthalmology, 90, 262– 267. Rouhi, R., Jafari, M., Kasaei, S., & Keshavarzian, P. (2015). Benign and malignant breast tumors classification based on region growing and cnn segmentation. Expert Systems with Applications, 42, 990 – 1002. dos Santos Ferreira, M. V., de Carvalho Filho, A. O., de Sousa, A. D., Silva, A. C., & Gattass, M. (2018). Convolutional neural network and texture descriptor-based automatic detection and diagnosis of glaucoma. Expert Systems with Applications, 110, 250 – 263. Shabut, A. M., Tania, M. H., Lwin, K. T., Evans, B. A., Yusof, N. A., AbuHassan, K. J., & Hossain, M. (2018). An intelligent mobile-enabled expert system for tuberculosis disease diagnosis in real time. Expert Systems with Applications, 114, 65 – 77. Soomro, T. A., Gao, J., Khan, T., Hani, A. F. M., Khan, M. A., & Paul, M. (2017a). Computerised approaches for the detection of diabetic retinopathy using retinal fundus images: a survey. Pattern Analysis and Applications, 20, 927–961. Soomro, T. A., Khan, M. A., Gao, J., Khan, T. M., & Paul, M. (2017b). Contrast normalization steps for increased sensitivity of a retinal image segmentation method. Signal, Image and Video Processing, 11, 1509–1517. Soomro, T. A., Khan, T. M., Khan, M. A., Gao, J., Paul, M., & Zheng, L. (2018). Impact of ica-based image enhancement technique on retinal blood vessels segmentation. IEEE Access, 6, 3524–3538. Soudani, A., & Barhoumi, W. (2019). An image-based segmentation recommender using crowdsourcing and transfer learning for skin lesion extraction. Expert Systems with Applications, 118, 400 – 410. Strobl, C., Boulesteix, A.-L., Zeileis, A., & Hothorn, T. (2007). Bias in random forest variable importance measures: Illustrations, sources and a solution. BMC bioinformatics, 8, 25. Walter, T., Klein, J.-C., Massin, P., & Erginay, A. (2002). A contribution of image processing to the diagnosis of diabetic retinopathy-detection of exudates in color fundus images of the human retina. IEEE transactions on medical imaging, 21, 1236–1243. Welfer, D., Scharcanski, J., Kitamura, C. M., Dal Pizzol, M. M., Ludwig, L. W., & Marinho, D. R. (2010). Segmentation of the optic disk in color eye fundus images using an adaptive morphological approach. Computers in Biology and Medicine, 40, 124–137. Yin, F., Liu, J., Ong, S. H., Sun, Y., Wong, D. W., Tan, N. M., Cheung, C., Baskaran, M., Aung, T., & Wong, T. Y. (2011). Model-based optic nerve head segmentation on retinal fundus images. In Engineering in Medicine and Biology Society, EMBC, 2011 Annual International Conference of the IEEE (pp. 2626–2629). IEEE. Yin, F., Liu, J., Wong, D. W. K., Tan, N. M., Cheung, C., Baskaran, M., Aung, T., & Wong, T. Y. (2012). Automated segmentation of optic disc and optic cup in fundus images for glaucoma diagnosis. In Computer-based medical systems (CBMS), 2012 25th international symposium on (pp. 1–6). IEEE. Zahoor, M. N., & Fraz, M. M. (2017). Fast optic disc segmentation in retina using polar transform. IEEE Access, 5, 12293–12300. Zhang, Z., Liu, J., Cherian, N. S., Sun, Y., Lim, J. H., Wong, W. K., Tan, N. M., Lu, S., Li, H., & Wong, T. Y. (2009). Convex hull based neuro-retinal optic cup ellipse optimization in glaucoma diagnosis. In Engineering in Medicine and Biology Society, 2009. EMBC 2009. Annual International Conference of the IEEE (pp. 1441–1444). IEEE. Zou, K. H., Warfield, S. K., Bharatha, A., Tempany, C. M., Kaus, M. R., Haker, S. J., Wells, W. M., Jolesz, F. A., & Kikinis, R. (2004). Statistical validation of image segmentation quality based on a spatial overlap index1: scientific reports. Academic radiology, 11, 178–189.

AN US

779

M

778

ED

777

PT

776

CE

775

Achanta, R., Hemami, S., Estrada, F., & Susstrunk, S. (2009). Frequency-tuned844 salient region detection. In 2009 IEEE Conference on Computer Vision and845 Pattern Recognition (pp. 1597–1604). 846 Achanta, R., Shaji, A., Smith, K., Lucchi, A., Fua, P., & S¨usstrunk, S. (2012).847 Slic superpixels compared to state-of-the-art superpixel methods. IEEE848 transactions on pattern analysis and machine intelligence, 34, 2274–2282. 849 Almazroa, A., Burman, R., Raahemifar, K., & Lakshminarayanan, V. (2015).850 Optic disc and optic cup segmentation methodologies for glaucoma image851 detection: a survey. Journal of ophthalmology, 2015. 852 Bernal, J., Kushibar, K., Asfaw, D. S., Valverde, S., Oliver, A., Mart´ı, R.,853 & Llad´o, X. (2017). Deep convolutional neural networks for brain im-854 age analysis on magnetic resonance imaging: a review. arXiv preprint855 arXiv:1712.03747, . 856 Boser, B. E., Guyon, I. M., & Vapnik, V. N. (1992). A training algorithm for857 optimal margin classifiers. In Proceedings of the fifth annual workshop on858 Computational learning theory (pp. 144–152). 859 Breiman, L. (2001). Random forests. Machine learning, 45, 5–32. 860 Canny, J. (1987). A computational approach to edge detection. In Readings in861 Computer Vision (pp. 184–203). Elsevier. 862 Cao, X. H., Stojkovic, I., & Obradovic, Z. (2016). A robust data scaling algo-863 rithm to improve classification accuracies in biomedical data. BMC bioin-864 formatics, 17, 359. 865 Cheng, J., Liu, J., Wong, D. W. K., Yin, F., Cheung, C., Baskaran, M., Aung,866 T., & Wong, T. Y. (2011). Automatic optic disc segmentation with peripap-867 illary atrophy elimination. In Engineering in Medicine and Biology Soci-868 ety, EMBC, 2011 Annual International Conference of the IEEE (pp. 6224–869 6227). IEEE. 870 Fan, Z., Rong, Y., Cai, X., Lu, J., Li, W., Lin, H., & Chen, X. (2018). Optic871 disk detection in fundus image based on structured learning. IEEE journal872 of biomedical and health informatics, 22, 224–234. 873 Freund, Y., & Schapire, R. E. (1997). A decision-theoretic generalization of874 on-line learning and an application to boosting. Journal of computer and875 system sciences, 55, 119–139. 876 Fu, H., Cheng, J., Xu, Y., Wong, D. W. K., Liu, J., & Cao, X. (2018). Joint877 optic disc and cup segmentation based on multi-label deep network and polar878 transformation. arXiv preprint arXiv:1801.00926, . 879 Gonzalez, R. (2006). Re woods, digital image processing, dorling kindersley. 880 Harizman, N., Oliveira, C., Chiang, A., Tello, C., Marmor, M., Ritch, R., &881 Liebmann, J. M. (2006). The isnt rule and differentiation of normal from882 883 glaucomatous eyes. Archives of ophthalmology, 124, 1579–1583. Henriksen, J. J. (2007). 3d surface tracking and approximation using gabor884 filters. South Denmark University, 28. 885 Jain, A. K. (1989). Fundamentals of digital image processing (prentice hall886 information and system sciences series). 887 Jonas, J. B., Fern´andez, M. C., & Naumann, G. O. (1992). Glaucomatous para-888 papillary atrophy: occurrence and correlations. Archives of Ophthalmology,889 110, 214–222. 890 Khan, M. A., Khan, T. M., Bailey, D., & Soomro, T. A. (2018). A general-891 ized multi-scale line-detection method to boost retinal vessel segmentation892 sensitivity. Pattern Analysis and Applications, (pp. 1–20). 893 Khan, M. A., Khan, T. M., Soomro, T. A., Mir, N., & Gao, J. (2017). Boosting894 sensitivity of a retinal vessel segmentation algorithm. Pattern Analysis and895 Applications, (pp. 1–17). 896 Lalonde, M., Beaulieu, M., & Gagnon, L. (2001). Fast and robust optic disc de-897 tection using pyramidal decomposition and hausdorff-based template match-898 ing. IEEE transactions on medical imaging, 20, 1193–1200. 899 Lee, R., Wong, T. Y., & Sabanayagam, C. (2015). Epidemiology of diabetic900 retinopathy, diabetic macular edema and related vision loss. Eye and vision,901 2, 17. 902 Liao, P.-S., Chen, T.-S., Chung, P.-C. et al. (2001). A fast algorithm for multi-903 level thresholding. J. Inf. Sci. Eng., 17, 713–727. 904 Lowell, J., Hunter, A., Steel, D., Basu, A., Ryder, R., Fletcher, E., & Kennedy,905 L. (2004). Optic nerve head segmentation. IEEE Transactions on medical906 Imaging, 23, 256–264. 907 Maldhure, M. P. N., & Dixit, V. (2015). Glaucoma detection using optic cup908 and optic disc segmentation. International Journal of Engineering Trends909 and Technology (IJETT), 20, 52–55. 910 Maninis, K., Pont-Tuset, J., Arbel´aez, P., & Gool, L. V. (2016a). Deep retinal911 image understanding. In Medical Image Computing and Computer-Assisted912 Intervention (MICCAI). 913 Maninis, K.-K., Pont-Tuset, J., Arbel´aez, P., & Van Gool, L. (2016b). Deep

AC

773 774

16

ACCEPTED MANUSCRIPT

914 915 916 917

Author Contribution Statement Zaka Ur Rehman: Conceptualization; Data curation; Formal analysis; Investigation; Methodology; Validation; Visualization; Roles/Writing - original draft; Writing - review & editing.

918 919 920 921

Supervision; Project administration; Resources; Software; Funding acquisition; Syed S. Naqvi: Supervision; Project administration; Resources; Writing - review & editing.

922

924 925 926 927

Tariq M. Khan: Supervision; Project administration; Resources; Writing - review & editing. Muhammad Arsalan: Visualization, Writing - review & editing Muhammad A. Khan: Writing - review & editing

CR IP T

923

928

CE

PT

ED

M

AN US

M. A. Khalil: Writing - review & editing

AC

929

17