Abstract

Glaucoma is a group of eye conditions, which can seriously damage optic nerves because of an elevated intraocular pressure. Nowadays, glaucoma has become one of the principal causes of blindness that results in irreversible visual loss. Early screening and treatment of glaucoma can prevent further progression of optic nerve degeneration effectively. The vertical cup-to-disc ratio (CDR) is a commonly used measurement for the detection of glaucoma, and therefore accurate segmentation of optic disc (OD) and optic cup (OC) regions in retinal fundus images is of great significance. In this paper, we present a prior shape constraint-based two-layer level set method for OD and OC segmentation in fundus images. This method uses two different layers of one level set function to represent the OD and OC boundaries. In this method, the distance regularized term is designed to guarantee that the distance between the OD and OC varies smoothly. By introducing the prior shape constraint term energy, the final segmentation results of OD and OC are always in the shape of approximate ellipses. In addition, the proposed method has the property of dealing with the intensity inhomogeneity of fundus images through the local fitting energy embedded. Experiments on images from the Baidu Research database demonstrate that the proposed method has superior performance in terms of accuracy and effectiveness.

1. Introduction

Nowadays, eye diseases are serious diseases that damage human health. The most common eye diseases include glaucoma, diabetic retinopathy, and maculopathy [1]. Because of causing optic nerve injury, the visual impairment caused by glaucoma is irreversible. Glaucoma becomes the second principal leading cause of blindness. Early detection and treatment of glaucoma can help prevent the vision from further loss and maintain the available visual function for most patients [2ā€“4]. In addition, the early screening of glaucoma is usually based on the length ratio of the vertical cup diameter to the vertical disc diameter in a retinal fundus image, which is widely used in clinical diagnosis of various eye diseases [5]. Figure 1 shows two fundus images with different cases, one is the normal case and the other is the glaucoma case. The optic disc (OD) is formed by the aggregation of retinal nerve fibers, which is light yellow in color and ellipse in shape. The optic cup (OC) is a cup-like area in the center of OD with variable sizes. The vertical cup-to-disc ratio (CDR) plays an important role in the diagnosis and classification of glaucoma [6, 7]. The larger the value of the CDR is, the greater the risk of glaucoma will be [8]. The effective CDR depends on the precise extraction of OD and OC regions in the fundus images. Therefore, accuracy and fast segmentation of OD and OC regions have attracted wide attention from researchers [6]. The traditional manual segmentation of OD and OC regions is of low efficiency and consumes a lot of manual effort [9]. In recent years, with the development of image processing and computer vision technology, automatic and accurate segmentation of OD and OC regions by computer-aided detection technology has become an important project [10].

Researchers have proposed a great number of methods for the segmentation of OD and OC. Generally, these methods fall into two broad categories: nonmodel-based methods and model-based methods [11]. The former ones mainly rely on image intensity, color, and other characteristics to segment the OD and OC, including the thresholding method, the pixel classification method, and the morphological processing method. The latter ones construct the constrain conditions to extract the OD and OC contours, which can be summarized as active shape model methods [12] and active contour model methods [13].

There are a large number of nonmodel algorithms proposed for the segmentation of OD and OC [11, 13]. In most cases, the boundary of the OD and OC are approximated as circular or elliptical shape. In [14], morphological operation and automatic thresholding methods were integrated for the segmentation of OD in fundus images. Firstly, in order to reduce the influence of blood vessels, a series of morphological open and close operations were applied to the original image to obtain a bright region-enhanced image. Secondly, a reduced region of interest was obtained by the automatic threshold method. Finally, the Hough transform was used to realize the segmentation of OD. In [15], the authors proposed an automatic OD localization and detection method. The first step was to locate the OD and obtain the region of interest (ROI). The second step was to detect the OD boundary by Canny operator in the ROI. In the third step, the threshold method was applied to the Hough transform image and the candidate circles were obtained. Finally, a comparison among candidate circles was made to extract the OD boundary. In [16], the K-means clustering method was applied to the segmentation of OD and OC. Firstly, the region of interest near the OD was extracted using an intensity-weighted centroid method. Secondly, through morphological operations in red channel and color space transform of the preprocessed RGB image, the color-enhanced image was obtained, on which the K-means clustering method was implemented to extract the OD areas. The similar method was used in OC segmentation, except that the preprocessing was performed in green channel. In [17], the fuzzy c-means method which was integrated with morphological techniques was applied to the segmentation of OD and OC and provided desirable results. In addition, all these operations were carried out in green channel of fundus images. Cheng et al. [18] presented a novel OD and OC segmentation method based on superpixel classification. The superpixels were determined to OD region or background by a contrast-enhanced histogram and a center surround statistics method, on the basis of which the OC region was estimated through the prior information of location.

Compared with the nonmodel-based methods, the model-based methods which involve active shape model methods and active contour model methods have better segmentation performance in terms of robustness and accuracy. The active shape model (ASM) describes an object in the original image by a statistical model of shape and appearance, which is used to identify the shape of an object in a new image [19]. In [20], Li and Chutatape proposed a modified ASM method for the extraction of OD boundary. The modified ASM method had better properties with robustness to the weak edge of the OD and vascular interference than the original ASM. In [21], Yin et al. proposed a method in which the ASM, edge detection, and circular Hough transform were utilized together to detect the OD in fundus images. The edge detection and circular Hough transform were applied to detect the edge of the OD as the initial contour of the ASM method. In [22], Cheng et al. proposed a novel ASM method based on sparse dissimilarity constrained coding for OD segmentation and reconstruction and obtained desirable results. However, the ASM methods cannot always achieve satisfactory results in practical applications because the number of parameters involved in the training set is limited. Unlike the ASM methods, active contour model methods can be characterized by the integration of image data, initial estimation, and basic knowledge constraints into one feature extraction process. They have been extensively used in the segmentation of OD and OC in fundus images. In [23], Simon and Michael applied the snake model to detect the OD in fundus images for the first time. In [24], Zhou et al. improved the traditional gradient vector flow (GVF-snake) to achieve accurate segmentation of OD boundary. The mean shift operator was introduced to the standard GVF cost function which can drive the internal and external energies towards the correct direction. This method obtained desirable segmentation results. In [25], Wong et al. proposed an edge-based variational level set method with the ellipse fitting energy to extract the OD boundary in fundus images. Then, the OC region was obtained by a threshold level set method. In [26], Dai et al. proposed a novel active contour model that consists of boundary energy, shape energy, and region energy to segment the OD region automatically in fundus images. This method was robust to the noises at boundary and initial contours due to the multiple energies involved. In [27], Gao et al. introduced the priori elliptical shape constraint of the OD into the energy model, which can improve the accuracy of the segmentation algorithm. In [28], Mittapalli and Kande proposed a modified local binary energy fitting region active contour model to extract the OD contour, on the basis of which the OC region was segmented according to the structure and gray level properties of the optic cup.

As mentioned above, most of the active contour model methods can be considered as two-step algorithms, in which the segmentation of OD and OC are achieved by step. These methods cannot obtain the segmentation results of OD and OC simultaneously. The structural relationship between optic OD and OC in the fundus images is not fully considered. In addition, the segmentation algorithm is affected by several factors, such as retinal pathological regions, vascular coverage, and inhomogeneity of fundus images, which may cause undesirable results. In order to solve the problem of segmentation of OD and OC in fundus images, this paper proposes a distance regularized two-layer level set method based on prior elliptic shape constraint (DRLSEC). This method represents the OD and OC boundaries in a fundus image by two different layers of one level set function, which has the following three advantages. (1) The distance relationship between OD and OC is taken into account, and only one level set function is used to realize the segmentation of OD and OC regions simultaneously. (2) This method is able to deal with images with intensity inhomogeneity. (3) This method overcomes the influence of vascular structure on the segmentation results of OD and OC effectively.

The remaining paper is organized as follows. In Section 2, the framework of the DRLSEC algorithm for OD and OC segmentation is described. Section 3 is the numerical calculation of the proposed DRLSEC algorithm. In Section 4, experimental results and comparison experiments are presented. Finally, this paper is concluded in Section 5.

2. Methodology

2.1. Distance Regularized Two-Layer Level Set Segmentation Method

Feng et al. [29] proposed a distance regularized two-layer level set (DR2LS) method to extract the inside and outside contours of some circular objects using only one level set function. In this method, the inside and outside contours of the object are represented by the 0-level contour and the k-level contour of the level set function. Generally, the object exterior region, the circular objects, and the object interior region have different intensities. In some traditional methods, two level set functions are employed to segment such a circular object. Hence, the iterative solution of two level set evolution functions is needed to be carried out simultaneously to obtain the final segmentation result. This process is computationally expensive [30]. In contrast, this framework is capable of achieving much more computation efficient segmentation with two different levels representing object boundary.

Let be the image domain, and image is defined on the domain. is the level set function. The 0-level of the level set function is denoted by , i.e., , and the k-level of the level set function is denoted by , i.e., as shown in Figure 2. In domain , the image is divided into three regions by and , and these regions can be defined as follows:

The Heaviside function H is introduced to obtain the following representation of the regions , , and :where in the case of and in the case of .

In [29], the following energy functional was proposed, and is expressed aswhere is the data term, which can be designed according to the properties of images. is the distance regularization term, which can constrain the varying of distance between the 0-level contour and the k-level contour of the level set function. The length regularization term is applied to derive the smooth 0-level and k-level contours by penalizing their length.

The data term is usually flexible. The data term from some traditional active contour models all can be used in this framework. Feng et al. [29] suggested the following data term from the famous region scalable fitting (RSF) model [31] to deal with intensity inhomogeneity of images, which is given bywhere is the function which fits values of local intensities and is the positive constant as weight coefficient. The function is a truncated kernel function with the following form:where to facilitate .

In order to make the distance between 0-level contour and k-level contour of the level set function vary slowly and smoothly, it is necessary to establish a constraint relationship between them. The distance regularization term D is defined as follows:where and are weighting positive constants. In the first term of the formula, is forced to be a smooth function . And the second term is employed to guarantee the function to be smooth.

The length regularization term is defined in the following form:where and are weight coefficients. The length of the 0-level contour and k-level contour is penalized by the first term and the second term, respectively, by which the contours of the 0-level and k-level guarantee to be smooth during the whole evolution.

2.2. The Proposed Method

The DR2LS method is suitable for the extraction of object contour with circular structure. In [29], this method is used for contour extraction of endocardium and epicardium of the left ventricle. The contours of the OD and OC can be seen as a circular structure in fundus images. But unlike the characteristics of contours of endocardium and epicardium of the left ventricle, there are many lesion areas in OD region. In addition, a large amount of vascular structure exists near the boundary of the OC. Both of them interfere with contours extraction greatly. As a result, the OD and OC contours cannot be extracted by using the DR2LS method directly. Therefore, in this paper, we improve the DR2LS method for meeting the OC and OD segmentation requirements.

Generally, the fundus images are in the presence of intensity inhomogeneities. Though the DR2LS method has suggested a data term to deal with intensity inhomogeneities, such a data term which only contains the local image information can be trapped into local minima easily and fails to provide a desirable result. Therefore, in this paper a more robust data term is strongly required to overcome this problem.

More specifically, an observed image can be formulated as follows:where denotes real true image which is considered as piecewise constant, denotes the bias field and varies slowly, and is the noise of normal distribution with zero mean. Based on this model, Li et al. [32] improved the data term of the RSF model with the following data term:

In the improved data term, replaces which means both the global and local image information are utilized to fit the local image intensities. This improvement can effectively avoid the problem associated with being tapped into local minima. Therefore, we employed it as the data term in the proposed method.

The OD has a slightly vertically elliptical shape and is divided into two separate regions: the central OC region with horizontally ellipse shape and neuroretinal rim [9, 33]. There are a large number of lesion areas and vascular structures in fundus image, which may cover some parts of the OC and OD contours. In order to overcome this problem, we introduced the prior elliptic shape constraint of the OD and OC into the framework of the DR2LS method, which can always force the contours of 0-level and k-level of the level set function to approximate to ellipse during the evolution. Such constraints can eliminate the influence of vascular structures and lesion areas on contours extraction.

The prior elliptic shape constraint of optic cup contour is defined as follows:and the prior elliptic shape constraint of optic disc contour is defined as follows:where and are the central coordinates of the two ellipses; and are the semimajor axis length; and are the semiminor axis length; and and are angles of rotation. For the convenience of discussion, the parameters involved in the prior elliptic shape constraint are represented as (, , , , ) in this paper. and are level set functions for describing the ellipse shape. The prior shape constraint can be defined as follows:where is the Heaviside function mentioned above and and are the weight coefficients. The first term is the shape constraint of OC, and the second term is the shape constraint of OD.

We introduce the new data term and the prior shape constraint term into the distance regularized two-layer level set segmentation framework, which is named DRLSEC, and its integral energy functional is in the following form:

The extraction of OD and OC contours can be achieved by optimizing the level set function such that the 0-level and k-level of the level set function can extract the OD and OC contours accurately.

3. Numerical Calculation

We can minimize the energy functional E alternately. For the parameters , , , , , , , , , and fixed, we use the standard gradient descent method to minimize the energy functional with respect to and obtain the corresponding gradient descent flow equation:where is representative of , which can be rewritten as . is the derivative of the Heaviside function. is the divergence operator.

For the parameters , , , , , , , and fixed, the energy functional is minimized with respect to , , and , which are given as follows:

For the parameter fixed, the distance regularization term is minimized by using the gradient descent flow method. Therefore, the evolution equation with respect to can be obtained:

For the parameters , , , , , , , , , and fixed, we can compute the bias field by minimizing the energy functional , which can expressed as

For the parameters , , , , , and fixed, the evolution function of , , , , and can be obtained by minimizing the energy functional , which is given bywhere , , , and are in the following forms:

In practical implementation, the regularized Heaviside function is used to maintain the stability of numerical calculations, which is given byand its derivative is derived in the form below:

In this paper, we use an explicit scheme to solve the numerical problem. All the spatial partial derivatives and and are approximated by the central difference, and the temporal partial derivative is approximated by the forward difference [34]. Then, equation (13) can be discretized aswhere is the time step and is the number of iterations. The final segmentation result is obtained though the iterative solution of the formula above.

In the evolution process of the level set function, the reaction-diffusion method [35] is introduced to solve the reinitialization problem:

The specific steps of the proposed method are presented as follows:ā€‰Step 1: initialize level set functions , , and as signed distance functions; initialize the function and bias field .ā€‰Step 2: calculate , , and according to equations (15)ā€“(17).ā€‰Step 3: update , , , , and according to equations (20)ā€“(29).ā€‰Step 4: evolve the level set function according to equation (14).ā€‰Step 5: smooth the level set function using equation (34).ā€‰Step 6: calculate using equation (18).ā€‰Step 7: calculate using equation (19).ā€‰Step 8: verify the convergence condition of level set function. The following criterion is employed to measure whether the level set function converges or not, which is defined byā€‰where is the sum of length changes of 0-level contour and k-level contour between and iterations. When ( is a threshold), the level set function will stop the evolution. Otherwise, the algorithm returns to Step 2 and continues.

4. Experimental Results

The experimental results are obtained using MATLAB R2018a on a computer with Intel (R) Core (TM) i7 and 64-bit Windows 10. It is applied to segment the OD and OC regions in fundus images on the database from Baidu Research with the webpage (http://refuge.grand-challenge.org), which collects data from the various business units and provides these datasets at no cost for research and personal uses [36]. There are 800 retinal fundus images in the database with annotation for OD and OC given. We apply the proposed DRLSEC method, the multiphase C-V method [37], the distance regularized two-layer level set method (DR2LS) [29], and the proposed model without prior elliptic shape to carry out the segmentation of OD and OC regions. The comparisons between these experimental results are made to evaluate the segmentation performance of these methods. The experimental parameters are set as follows: , , , , , , , , , , , , , and . The Dice coefficient which represents the average overlapping ratio [38] is employed as a metric to evaluate the segmentation accuracy, which computes the similarity between two sets. The similarity between the object region Y obtained by a segmentation algorithm and the ground truth region X is computed. The Dice coefficient is given as

It can be easily observed that the value of the Dice coefficient ranges from 0 to 1, and the larger the Dice coefficient is, the higher is the accuracy of the algorithm.

In this section, we analyze the effect of the prior elliptic shape constraint item on the segmentation results. The comparison results are illustrated in Figure 3. The proposed method without prior elliptic shape constraint term and the proposed method are used to segment the OD and OC regions in fundus images. The experimental results are shown in 3(a) and 3(b), respectively. In each subgraph, the left column displays both the 0-level and k-level contours of the level set function. The right column displays the corresponding level set function in three-dimensional manner. In addition, we adopt the negative version of the level set function for visual display. The 0-level contour is represented by a curve in green and the k-level contour is in red. It can be observed from the evolution procedure of level set functions that the prior elliptic shape constraint term always forces both the 0-level and k-level contours to approximate ellipse shapes. Meanwhile, the level set functions are always smooth without reinitialization during the entire evolution process, which guarantee a stable numerical calculation for the accurate segmentation results.

We evaluate the performance of the proposed method in different stages of glaucomatous cupping. There are mild, moderate, and severe classes representing different stages of glaucoma. The proposed DRLSEC method is compared with the multiphase C-V method and distanced regularized two-layer level set (DR2LS) method, and the corresponding segmentation results are displayed in Figures 4ā€“6. In these figures, it is obvious that the proposed DRLSEC method and DR2LS method are able to deal with intensity inhomogeneity effectively, which is due to their local intensity fitting energy. In contrast, the multiphase C-V method which lacks local image information cannot overcome the problem associated with the intensity inhomogeneity of retinal fundus images, and it obtains undesirable segmentation results. Compared with the DR2LS method, the proposed DRLSEC method with the ellipse prior shape constraints can avoid the interference of vascular structures and lesion areas, and it can always gain the desirable OD and OC contours. As can be seen from the three groups experimental in different glaucoma stages, our method can overcome the intensity inhomogeneity and interference of vascular structures and lesion areas, obtaining desirable results superior to that of the other two methods.

In order to quantitatively evaluate the performance of the proposed method, we use the Dice coefficient as a criterion to estimate the accuracy of the segmentation results. We compute the average overlapping ratio of OD and OC according to the segmentation results of 800 images in the database. The proposed DRLSEC method is compared with the multiphase C-V method, DR2LS method, DR2LS method with EC (elliptic shape constraint), and DRLSEC method without EC, and the comparison results of all the methods are listed in Table 1. As shown in Table 1, the performance of the DRLSEC method is better than the other methods. It has obtained the average overlapping ratio of 72.16% and 66.45% in OD and OC segmentation, respectively.

5. Conclusions

In this paper, we have presented a distance regularized two-layer level set method based on prior elliptic shape constraint for the segmentation of the OD and OC regions in retinal fundus images. Accurate segmentation of OD and OC regions is conducive to achievement of the effective CDR which determines the diagnosis and classification of glaucoma. The proposed method mainly consists of three energy constraints. The improved local energy fitting term makes the algorithm have the intrinsic property of processing the intensity inhomogeneity of the fundus images. The distance regularization term constrains the varying of distance between the OD and OC. The elliptic shape constraints term is employed to force the final OD and OC segmentation results approximate to two ellipses. Experimental results and quantitative analysis have shown that the proposed method is able to segment the OD and OC regions simultaneously and outperforms the other three methods in robustness and accuracy.

Data Availability

The fundus image data used to support the findings of this study are available at http://refuge.grand-challenge.org.

Conflicts of Interest

The authors declare no conflicts of interest.

Acknowledgments

This research was supported in part by the National Natural Science Foundation of China under Grant nos. 61701101, U1713216, 61901098, and 61971118, the National Key Robot Project under Grant no. 2017YFB1301103, and the Fundamental Research Fund for the Central Universities of China N172603001, N181602014, N172604004, N172604003, and N172604002.