Automated Optic Disc Segmentation Using Basis Splines-Based Active Contour

We develop a fully automated method for the segmentation of optic disc in retinal fundus images using basis-spline-based active contour. The segmentation is achieved by performing scaling, translation, and rotation of the active contour, thereby giving rise to five free parameters. The energy of the active contour is defined by the local contrast and is optimized with respect to five free parameters to get the best fit on the optic disc using gradient descent technique and Green’s theorem. The use of gradient descent technique and Green’s theorem reduces the computational cost and speeds up the segmentation task. The detection of optic disc is achieved using multiresolution-based normalized cross-correlation technique. The detection point is used for the initialization of the active contour. Subsequent optimized evolution of the basis-spline-based active contour provides an accurate segmentation of the optic disc. We present validations on the databases such as Drishti-GS, MESSIDOR, RIGA, and a local database containing 101, 1200, 750, and 942 retinal fundus images, respectively, amounting to a total of 2993 fundus images. Their corresponding Dice index scores are 0.9182, 0.8912, 0.9331, and 0.9343. Basic data exploration is done on the results obtained to visualize the trends and distribution of the performance parameters throughout the databases. This also helps us evaluate the algorithm’s overall performance more accurately.


I. INTRODUCTION
Glaucoma is caused due to increased intraocular pressure inside the eye. If left untreated it may lead to partial or complete loss of vision. It is relatively common among adults over 60 years of age and is the second leading cause of blindness in the world. The loss of vision due to glaucoma cannot be recovered but can be slowed or prevented if recognized early [1]. It is estimated that over 3% of the global population over the age of 40 years is suffering from glaucoma, the majority are undiagnosed. In 2013, 64.3 million people, aged 40-80 years, were estimated to be suffering from glaucoma worldwide [1]. As per the World Report on Vision published by World Health Organization (WHO) in 2019, about The associate editor coordinating the review of this manuscript and approving it for publication was Vishal Srivastava.
2.2 billion people suffer from visual impairments. Among the 2.2 billion sufferings, at least 1 billion people could have been saved from blindness if they had received appropriate treatment on time [2]. The distribution of eye impairments is not even in the world. But on average, about 53.40% and 51.65% of people from Asia-Pacific and East Asian highincome regions suffer from myopia. As of 2020, a staggering population of 195.6 million suffer from age-related macular degeneration and a population of 76 million suffer from glaucoma. The number is predicted to rise to about 243.4 million and 95.4 million respectively, in the next decade [2].
As glaucoma progresses, the various geometrical parameters of the retina undergo visible changes. It is key to monitor the optic disc (OD) for the early detection and assessment of glaucoma. The OD is also called the optic nerve head. It is the circular region in the retina where the optic nerves meet. To the human eye, it can be approximated as being circular with its color ranging from orange to pink and white as well, as shown in Fig.1. Presently, ophthalmologists perform OD detection and outlining manually. This is quite an arduous task and the results for the same may vary from one ophthalmologist to another.
As an aid to the traditional OD segmentation methods and also as an objective way of segmenting OD, we propose an automated and accurate OD segmentation method.

A. OUR CONTRIBUTION
We propose a B-spline-based active contour model for automated segmentation of OD in retinal fundus images. Initialization of the B-spline-based active contour is achieved using the multiresolution-based normalized cross-correlation (MNCC) technique. X-plane of the XYZ color scheme is used to obtain better OD detection accuracy. The OD of the fundus image is then cropped to reduce the computations of segmentation. This reduces the time taken for segmentation considerably. The coordinates of detected OD are used for the initialization of B-spline-based active contour. The active contour evolves from a shape-specific (circle) initialization to an amorphous shape for accurate segmentation of the OD. Gradient descent and Green's theorem are then used to minimize the local energy function with respect to five free parameters. An exploratory data analysis of the segmentation is done to get a clear picture of the performance of the proposed method through the database.

II. PRIOR RESEARCH
Over the years, various state-of-the-art techniques, unique in theoretical modeling and approach have been brought about. Morales et al. proposed a watershed transformation based OD detection method. They used circle-fitting technique for the segmentation of OD. They used principal component analysis (PCA) with other morphological operations as a preprocessing step to improve the accuracy of OD detection [3]. Zahoor et al. utilized morphological operations and circular Hough transform for preprocessing and OD localization, respectively. The OD was segmented using the implementation of polar transform [4]. Cheng et al. developed a superpixel-based OD segmentation algorithm. They used histogram analysis with center-surround statistics to categorize the OD [5]. Dashtbozorg et al. segmented the OD by introducing a sliding-band filter [6]. Sigut et al. devised a contrastbased circular approximation method to detect and segment the OD [7]. Joshi et al. developed a region-based active contour and traditional Chan-Vese model for the segmentation of OD. It was a modified version of the active contour and traditional Chan-Vese model [8] with the modification being the use of the local image information from the area around the region of interest [9]. Kumar et al. developed a circular active contour model for the segmentation of OD [10], [11]. Abdullah et al. proposed a technique for localization and segmentation of OD in retinal images using circular Hough transform and grow-cut algorithm [12]. Dey et al. proposed affine snakes in gradient vector field technique for the segmentation of OD [13], [14]. Almazroa et al. developed a level set-based algorithm for the segmentation of OD. The blood vessel occlusions that interfere with the segmentation of OD were dealt with using an image inpainting technique [15]. Ramani et al. devised a pixel density-based OD detection method. They used circular Hough transform coupled with Hough peak value selection and red channel-based superpixel segmentation for OD segmentation [16].
Maninis et al. presented a Deep Retinal Image Understanding (DRIU) tool, which is a unified framework of retinal image analysis that provides both retinal vessel and OD segmentation [17]. Sevastopolsky devised a U-Net based OD segmentation algorithm [18]. It used the input image to generate a probability map. They used contrast limited adaptive histogram equalization (CLAHE) for preprocessing [19]. Zilly et al. developed a convolutional neural network (CNN) inspired ensemble learning technique for the segmentation of OD. Entropy sampling and Boosting were used to improve the developed network [20]. Mohan et al. proposed Fine-Net [21] and a combination of Fine-Net and P-Net [22] CNN models for the segmentation of OD. Wang et al. developed a coarse-to-fine U-Net model-based OD segmentation algorithm. They trained the network on a pair of color fundus images and their grayscale vessel density maps. They also attained two different segmentation results from the input image. This coarse segmentation obtained is then fed to a U-Net model for further fine segmentation [23]. Liu et al. devised an adversarial training method for the segmentation of OD. They used an improved U-Net-based algorithm for the detection of OD. The OD segmentation was achieved by a Patch-level adversarial network to enhance higher consistency between ground truth and algorithm's segmentation output [24]. Kadambi et al. developed a WGAN domain adaptation framework for the detection of OD and optic cup boundaries from the fundus images. The model is an adversarial domain adaptation framework guided by Wasserstein distance to improve stability and convergence [25]. Sun et al. have devised a deep object detection networkbased pipeline for the detection of OD. The OD boundary is then determined using the faster R-CNN method as VOLUME 10, 2022 an object detector by transforming the predicting bounding box into a vertical and non-rotated ellipse [26]. Zilly et al. proposed an ensemble learning-based CNN, where selective informative points are selected using an entropy sampling technique. A novel learning framework for the convolutional framework based on boosting is designed with the sampled points [27]. Zhang et al. developed a particle swarm optimization enhanced deep neural network algorithm for the segmentation of OD [28]. Fu et al. proposed a deep learning architecture, M-net, to segment OD and optic cup jointly as a one-stage multi-label system. Jiang et al. proposed an end-to-end region-based CNN for joint OD and optic cup segmentation. They also used Atrous convolution to boost the feature extraction performance [29]. Li et al. studied a sequence of supervised descent directions of the OD boundary coordinates. They also introduced histograms of gradient orientations to represent the OD. They felt that the performance of the OD segmentation can be improved by better shape-appearance modeling [30]. Thakur et al. approach. This method is used to segment the OD in abnormal fundus images. This model was proposed to be used when there is a lack of sufficient training samples. The probability bubble is dependent on the position relationship between retinal vessels of OD, and the density of the intersection points of the line segments that are fitted on the main blood vessels through Hough transforms [33]. Xiong et al. have introduced a Bayesian U-Net model for the segmentation of fundus images. They have attempted to address the limitations that arise due to pixel-level OD annotation mask, which induces inter-subject variance. They have addressed the limitations by introducing a weak label-based Bayesian U-Net that uses Hough transform-based annotations to segment OD in fundus images. The optimization of the proposed algorithm was achieved using Expectation-maximization algorithm. The authors claim that the performance of the proposed method is superior to that of fully and weakly supervised methods [34]. Pachade et al. developed a nested EfficientNet patch-based adversarial learning framework for the segmentation of OD and optic cup. The developed method uses the EfficientNetB4 as an encoder with a nested network of pre-activated residual blocks, atrous spatial pyramid pooling blocks, and attention gates. To obtain accurate segmentation the network is guided by a combination of cross-entropy and Dice coefficient loss [35]. Hasan et al. developed a DRNet model with a redesigned skip connection and called it residual skip connection. The developed DRNet model has fewer parameters, leading to shorter training and testing times. For localization of OD and fovea centres a 2D bell-shaped heatmap is used [36]. They used B-spline, a traditional image segmentation method, with statistical models to select the appropriate shape template for a more accurate segmentation [40]. This work is also motivated by the seminal contributions from Thevenaz et al. [41], Unser et al. [42], [43], Delgado-Gonzalo et al. [44], Mogali et al. [45], and Pediredla et al. [46].

III. B-SPLINE-BASED ACTIVE CONTOUR FORMULATION
We have used two concentric circles to construct the shape template model. The parameterization of the concentric shape template is done with the use of a cubic B-spline kernel. Its properties of non-negativity, compact support, continuous second-order derivative, and minimum curvature make it the prime choice [47], [48]. Additionally, it satisfies the Riesz basis property [49], which ensures the stability of representation between the continuous contour and the discrete coefficients used in the representation.

A. SPECIFICATION OF SHAPE TEMPLATE
Let r 0 (t) = (u 0 (t), v 0 (t)) T and r 1 (t) = (u 1 (t), v 1 (t)) T be a pair of parameterized inner and outer contours, respectively.  Here 't' is an independent variable. 0 and 1 are subscripts used to differentiate between the inner and outer contour, respectively. Kernel function ζ (t) and its integer shifted version is used to represent the function u 0 (t), u 1 (t), v 0 (t), and v 1 (t). The parameterization has been achieved in accordance with [50], and vector representation form of the results are: where and v i (t) are periodic for closed curves which leads to the following equivalent representation. where , are the periodized coefficient sequences with period M (M represents the number of knots). An example of a shape template designed using cubic B-spline with eight knots is shown in Fig. 2.
Outer contour (r 1 (t)) of an origin-centered convex shape can be a scaled variant of the inner contour (r 0 (t)), i.e., r 1 (t) = αr 0 (t), where α is a scalar variable that is greater than 1. As a result, the outer contour coefficients get scaled up by a factor of α times the inner contour coefficients. The local energy function's computation region and the control over the annular region width are achieved by the parameter α.

B. THE DYNAMIC ACTIVE CONTOUR
The derivation of active contours from the shape template in (2) is according to the following affine transformation: where i = 0, 1 represents inner and outer contours, respectively, of the active contour. For conciseness of notation, we have dropped the parameter t and denoted (U (t), V (t)) and (u(t), v(t)) as (U , V ) and (u, v), respectively. X 1 and X 2 are the scale parameters; θ is the angle of rotation; u c and v c are the translation parameters.

C. THE ENERGY FUNCTION
For an active contour initialized on an image p, let 0 and 1 be the regions enclosed by the inner and outer contours, respectively. We define the active contour energy as a normalized local contrast function: where a 0 and a 1 are the scalars (in the ratio 2:1) for the areas enclosed by the inner and outer contour and E 0 and E 1 represents the energy of the inner and outer contour of the shape template, respectively. This ensures the active contour is inert in regions of constant intensity. The term X 1 X 2 is the normalizing factor (net scaling in the area) of the active contour. Minimizing E enables the active contour to lock on to objects that are brighter than their immediate surroundings and maximizing it would have the opposite effect.

D. PARAMETERIZATION OF B-SPLINE COEFFICIENT ADAPTATION
Considering image p, finding partial derivatives of E with respect to the free parameters and expressing them in terms of u 0 (t), v 0 (t), u 1 (t), and v 1 (t), we have The weighted sum of B-splines, as in (1), can be used to express the active contour template. B-spline coefficient adaptation from the partial derivatives of parameters shown in (5) is done as in (6) and (7), shown at the bottom of the next page. Convergence towards the OD is accelerated due to the use of pre-computed kernels b k−l and ζ (t). The partial derivatives of the parameter are bound to be finite-length sequences as kernels b k−l and ζ (t) are finite.
The Q p,1 and Q p,2 are the common terms present in partial derivatives with respect to (X 1 , X 2 , and θ) and in partial derivatives with respect to (u c and v c ), respectively. The iterative computation of kernels Q p,1 and Q p,2 reduces the overall computations required. Recalculation of free parameters and change of the corresponding energy terms in the cost function is required after every iteration. The Locality property of B-splines means that its coefficients often remain unchanged.
Optimization of the partial derivatives of energy E with respect to the free parameters is achieved using gradient descent technique [51]. By optimizing energy E, the contrast between the two regions increases, in other words, the contrast in the annular region decreases.
The equations associated with gradient descent are where P n denotes the parameter (X 1 , X 2 , θ, u c , and v c ) at iteration 'n'. Green's theorem [52] is used to reduce the computations required to find the partial derivatives of the B-spline coefficient with respect to the free parameters.

E. AUTOMATIC LOCALIZATION
B-spline-based active contour is initialized on the retinal fundus image using the MNCC technique [53]. We calculate a normalized sliding inner-product between the scaled down version of the input image f , and the template h. The template is developed by cropping the OD region manually from a retinal fundus image.
To obtain better localization accuracy, X-plane of the XYZ color scheme is used. The comparison of Red channel from the RGB color scheme and X-plane of the XYZ color scheme is shown in Fig. 3. The X-plane image provides us a distinguishable and sharp OD area in red-dominated, well-lit and dull fundus images. The Red channel image fails to provide a shape and distinguishable OD region in red-dominated and dull fundus images. These factor influenced our decision to use the X-plane over the Red channel for the localization and segmentation of OD.
For fast localization, NCC is performed at the lowest level of a four-level subsampled pyramid [54] representation of the input fundus image and the OD location information is propagated to the highest level of the decomposition. The cosine similarity measure ξ between the input image f and the template h at the j th level of the pyramid defined by f (j) and h (j) is given by The location corresponding to the peak of the cosine similarity measure is chosen for initializing the B-spline-based active contour. An illustration of OD detection for initialization of the active contour is shown in Fig. 4.

F. ALGORITHM IN BRIEF
The OD segmentation using B-spline-based active contour algorithm, in brief, is given below: Input: Color fundus image Output: OD segmentation on the input color fundus image Initialization: Assign the path of the OD template image to be used by the MNCC algorithm. 1: Perform OD detection using MNCC algorithm 2: Crop the OD region from the fundus image using the detection coordinates obtained from Step 1. 3: Initialize B-spline-based active contour in the detection coordinates in the OD cropped image. 4: Perform affine transformation of the B-spline-based active contour minimizing the energy function. 5: Optimize the free parameters using gradient descent and Green's theorem for the best fit active contour on the OD. 6: Stitch the segmentation mask obtained using Step 5 to the original fundus image.

IV. EXPERIMENTS
The proposed method was validated on Drishti-GS [55], MESSIDOR [56], RIGA [57], and a locally acquired database totaling 2993 fundus images. Drishti-GS, MESSIDOR, and RIGA are publicly available datasets. All the databases consisted of colour fundus images. Drishti-GS database consists of 101 images captured from a population belonging to Indian ethnicity. Each image has 88156 VOLUME 10, 2022 a resolution of 2896 × 1944 pixels. The data collection and annotations were conducted by Aravind Eye Hospital, Madurai, India. This dataset is divided into two namely, the training and testing set consisting of 50 and 51 images respectively. Ground truth was collected from four experts with varying clinical experience. The database provides OD segmentation soft-maps fused on one binary image and average OD boundary derived from four expert markings.
MESSIDOR database is one of the most popular fundus imaging database consisting of 1200 images with resolutions of 2304 × 1536, 2240 × 1488, and 1440 × 960 pixels. The database provides OD ground truth and fovea center annotation by a single clinician. MESSIDOR stands for Methods to Evaluate Segmentation and Indexing Techniques in the field of Retinal Ophthalmology (in French). It was part of a research project conducted by the French Ministry of Research and Defense in the year 2004. The images were captured using a Topcon TRC NW6 non-mydriatic retinograph. Among the 1200 images, 800 were captured with pupil dilation and the rest were captured without pupil dilation.
RIGA is a de-identified database derived from three different sources namely, MESSIDOR, Bin Rushed eye center, and Magrabi eye center consisting of 460, 195, and where b k−l (t) = ζ (t)ζ (t + k − l) and M i denotes the number of knots. VOLUME 10, 2022 95 images respectively. The images were manually annotated by six experienced ophthalmologists. The database consists of annotations of OD and optic cup. The local database consists of 942 images captured from three different fundus imaging devices (Courtesy: Spectrum Lab, Department of Electrical Engineering, Indian Institute of Science, Bangalore, India). The database is provided with OD annotations by five expert ophthalmologists along with average OD ground truth.

V. RESULTS
We developed an ImageJ [58] plugin for the implementation of the proposed method. ImageJ is a popular Javabased biomedical image processing tool developed by W. Rasband and team at the National Institutes of Health (NIH), USA [59]. For validation, we have used the databases Drishti-GS [55], MESSIDOR [56], RIGA [57], and a local database amounting to a total of 2993 images. The segmentation of the OD using the proposed method on randomly selected fundus images from different databases have been shown in Fig. 5.
A visual representation to understand the segmentation performance metrics are shown in Fig. 6. Expert ground truth and algorithm segmentation are overlaid to determine the parameters such as true positives (green), true negatives (black), false positives (yellow), and false negatives (red). With most of the area representing the true positives (green), very negligible areas representing the false positives (yellow) and false negatives (red), it is evident that the expert ground truth and algorithm segmentation are in close agreement with each other.
Using the number of true positives, true negatives, false positives, and false negatives, we have evaluated the sensitivity (S e ), specificity (S p ), accuracy (A c ), and error (E r ). Our results are then compared against the state-of-the-art methods. We have compared the segmentation results with the ground truth and computed the Dice and Jaccard similarity indices. The Dice and Jaccard indices are evaluated using the following expressions: where D i and J i are the Dice and Jaccard similarity indices, respectively. X and Y are the segmentation result obtained from the proposed technique and the ground truth, respectively. The comparison of the Dice and Jaccard index variations within the databases is shown in Fig. 7. The comparison of the distribution of the Dice and Jaccard indices between databases is shown as box plots in Fig. 8. The comparison of the Dice index obtained for the six experts' annotations of the RIGA database are represented graphically in Fig. 9.
From the graph, we see that the Dice index value hovers around the 0.95 mark for all the annotations. So, we can infer that the proposed algorithm provides an accurate segmentation of the OD.  Table 1 summarizes the comparison of the performance parameters of the existing state-of-the-art techniques with the proposed algorithm.

VI. DISCUSSION
The proposed method is developed using traditional image processing concepts. The developed methods use the concepts such as multi-resolution-based normalized crosscorrelation method for the detection and basis splines-based active contour for the segmentation of the OD. Normalization and histogram equalization is used as a standard preprocessing step to improve the detection and segmentation. After detection of the OD but, prior to its segmentation, the image is cropped to retain only the region of interest to reduce computation and improve the execution speed.
Recently, various machine learning-based methods for the detection and segmentation of the OD have been reported. These methods have been included for comparison with the proposed method. From Table 1, it is evident that machine learning algorithms can perform better than traditional image processing methods when there is adequate training data available. But, it is difficult to acquire large amounts of medical imaging data with expert annotations. With the data currently available, machine learning algorithm's results appear to outperform traditional image processing methods marginally. Hence, it is our opinion that image processing methods are more robust when the available databases consist of a small number of images.
Although with data augmentation it might be possible for machine learning algorithms to outperform traditional image processing methods, there is a possibility it can lead to over-fitting. In which case the machine learning methods might perform poorly in a test set. These problems will not  be encountered in the case of traditional image processing methods.

VII. CONCLUSION
We proposed a novel basis spline-based active contour method for the segmentation of optic disc. The technique provides an accurate and amorphous segmentation of the optic disc against many of the techniques based on shape-specific segmentation. Hence, the segmentation outline provided by the proposed method is very close to the natural shape of the optic disc, which is amorphous. Optimization of energy function of basis spline-based active contour was achieved using gradient descent and Green's theorem. The developed algorithm was validated on three publicly available databases, Drishti-GS, MESSIDOR, RIGA, and one locally procured database amounting to a total of 2993 fundus images. The proposed technique achieves an overall sensitivity, specificity, accuracy, error rate, Jaccard, and Dice index scores of 94.07%, 99.82%, 99.71%, 0.28%, 85.59%, and 93.01% respectively. The performance metrics obtained prove that the proposed method has high accuracy in the highly varied fundus image data.