# Automatic Segmentation of Optic Disc and Cup in Retinal Fundus Images Using Improved Two-Layer Level Set Method

Authors: Ying Wang; Xiaosheng Yu; Jianning Chi; Chengdong Wu
Journal: Mathematical Problems in Engineering
 (2019)

---

## 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.


---

## Body

## 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].Figure 1
The area surrounded by the black curve is the optic cup area; the optic disc area is indicated by the red curve; and the neuroretinal rim is the area between them.Vdisc denotes the vertical optic disc diameter, and Vcup denotes the vertical optic cup diameter. (a) Normal case. (b) Glaucoma case.
(a)
(b)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. InSection 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Ω⊂R2 be the image domain, and image ϕ:Ω⟶R is defined on the domain. I:Ω⟶R is the level set function. The 0-level of the level set function ϕ is denoted by C0, i.e., C0≜x:ϕx=0, and the k-level of the level set function ϕ is denoted by Ck, i.e., Ck≜x:ϕx=0 as shown in Figure 2. In domain Ω, the image is divided into three regions by C0 and Ck, and these regions can be defined as follows:(1)Ω1≜x:ϕx<0,Ω2≜x:0<ϕx<k,Ω3≜x:ϕx>k.Figure 2
The area divided by the 0-level contour andk-level contour.The Heaviside functionH is introduced to obtain the following representation of the regions Ω1, Ω2, and Ω3:(2)M1ϕx=1−Hϕx,M2ϕx=Hϕx−Hϕx−k,M3ϕx=Hϕx−k,where Miϕx=1 in the case of x∈Ωi and Miϕx=0 in the case of x∉Ωi.In [29], the following energy functional Eseg was proposed, and is expressed as(3)Eseg=ε+D+L,where ε is the data term, which can be designed according to the properties of images. D 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 L 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 by(4)εϕ,f1,f2,f3=∑i=13λi∫∫Kρx−yIx−fiy2Miϕxdxdy,where fiyi=1,2,3 is the function which fits values of local intensities and λii=1,2,3 is the positive constant as weight coefficient. The function Kρ is a truncated kernel function with the following form:(5)Kρu=a,foru≤ρ,0,foru>ρ,where a>0 to facilitate ∫u≤ρKρu=1.In order to make the distance between 0-level contour andk-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:(6)Dϕ,α=μ1∫12∇ϕx−αx2dx+μ2∫12∇αx2dx,where μ1 and μ2 are weighting positive constants. In the first term of the formula, ∇ϕx is forced to be a smooth function αx. And the second term is employed to guarantee the function αx to be smooth.The length regularization termL is defined in the following form:(7)Lϕ=v1∫∇Hϕxdx+v2∫∇Hϕx−kdx,where v1>0 and v2>0 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:(8)Ix=bxcx+nx,where cx denotes real true image which is considered as piecewise constant, bx denotes the bias field and varies slowly, and nx 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:(9)εϕ,c1,c2,c3,bx=∫∑i=13λi∫Kρx−yIx−byci2Miϕxdxdy.In the improved data term,byci replaces fiy 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:(10)ϕ1x=1−x−x1cosθ1+y−y1sinθ12a12+−x−x1sinθ1+y−y1cosθ12b121/2,and the prior elliptic shape constraint of optic disc contour is defined as follows:(11)ϕ2x=1−x−x2cosθ2+y−y2sinθ22a22+−x−x2sinθ2+y−y2cosθ22b221/2,where x1,y1 and x2,y2 are the central coordinates of the two ellipses; a1 and a2 are the semimajor axis length; b1 and b2 are the semiminor axis length; and θ1 and θ2 are angles of rotation. For the convenience of discussion, the parameters involved in the prior elliptic shape constraint are represented as (a, b, θ, x, y) in this paper. ϕ1x and ϕ2x are level set functions for describing the ellipse shape. The prior shape constraint can be defined as follows:(12)Eellipse=γ1∫ΩHϕx−Hϕx2d+γ2∫ΩHϕx−k−Hϕx2dx,where Hx is the Heaviside function mentioned above and γ1>0 and γ2>0 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:(13)Eϕ,α,c1,c2,c3,bx,a,b,θ,x,y=Dϕ,α+εϕ,c1,c2,c3,bx+Lϕ+γ1∫ΩHϕx−Hϕ1x2dx+γ2∫ΩHϕx−k−Hϕ2x2dx.The extraction of OD and OC contours can be achieved by optimizing the level set function such that the 0-level andk-level of the level set function can extract the OD and OC contours accurately.


## 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Ω⊂R2 be the image domain, and image ϕ:Ω⟶R is defined on the domain. I:Ω⟶R is the level set function. The 0-level of the level set function ϕ is denoted by C0, i.e., C0≜x:ϕx=0, and the k-level of the level set function ϕ is denoted by Ck, i.e., Ck≜x:ϕx=0 as shown in Figure 2. In domain Ω, the image is divided into three regions by C0 and Ck, and these regions can be defined as follows:(1)Ω1≜x:ϕx<0,Ω2≜x:0<ϕx<k,Ω3≜x:ϕx>k.Figure 2
The area divided by the 0-level contour andk-level contour.The Heaviside functionH is introduced to obtain the following representation of the regions Ω1, Ω2, and Ω3:(2)M1ϕx=1−Hϕx,M2ϕx=Hϕx−Hϕx−k,M3ϕx=Hϕx−k,where Miϕx=1 in the case of x∈Ωi and Miϕx=0 in the case of x∉Ωi.In [29], the following energy functional Eseg was proposed, and is expressed as(3)Eseg=ε+D+L,where ε is the data term, which can be designed according to the properties of images. D 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 L 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 by(4)εϕ,f1,f2,f3=∑i=13λi∫∫Kρx−yIx−fiy2Miϕxdxdy,where fiyi=1,2,3 is the function which fits values of local intensities and λii=1,2,3 is the positive constant as weight coefficient. The function Kρ is a truncated kernel function with the following form:(5)Kρu=a,foru≤ρ,0,foru>ρ,where a>0 to facilitate ∫u≤ρKρu=1.In order to make the distance between 0-level contour andk-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:(6)Dϕ,α=μ1∫12∇ϕx−αx2dx+μ2∫12∇αx2dx,where μ1 and μ2 are weighting positive constants. In the first term of the formula, ∇ϕx is forced to be a smooth function αx. And the second term is employed to guarantee the function αx to be smooth.The length regularization termL is defined in the following form:(7)Lϕ=v1∫∇Hϕxdx+v2∫∇Hϕx−kdx,where v1>0 and v2>0 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:(8)Ix=bxcx+nx,where cx denotes real true image which is considered as piecewise constant, bx denotes the bias field and varies slowly, and nx 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:(9)εϕ,c1,c2,c3,bx=∫∑i=13λi∫Kρx−yIx−byci2Miϕxdxdy.In the improved data term,byci replaces fiy 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:(10)ϕ1x=1−x−x1cosθ1+y−y1sinθ12a12+−x−x1sinθ1+y−y1cosθ12b121/2,and the prior elliptic shape constraint of optic disc contour is defined as follows:(11)ϕ2x=1−x−x2cosθ2+y−y2sinθ22a22+−x−x2sinθ2+y−y2cosθ22b221/2,where x1,y1 and x2,y2 are the central coordinates of the two ellipses; a1 and a2 are the semimajor axis length; b1 and b2 are the semiminor axis length; and θ1 and θ2 are angles of rotation. For the convenience of discussion, the parameters involved in the prior elliptic shape constraint are represented as (a, b, θ, x, y) in this paper. ϕ1x and ϕ2x are level set functions for describing the ellipse shape. The prior shape constraint can be defined as follows:(12)Eellipse=γ1∫ΩHϕx−Hϕx2d+γ2∫ΩHϕx−k−Hϕx2dx,where Hx is the Heaviside function mentioned above and γ1>0 and γ2>0 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:(13)Eϕ,α,c1,c2,c3,bx,a,b,θ,x,y=Dϕ,α+εϕ,c1,c2,c3,bx+Lϕ+γ1∫ΩHϕx−Hϕ1x2dx+γ2∫ΩHϕx−k−Hϕ2x2dx.The extraction of OD and OC contours can be achieved by optimizing the level set function such that the 0-level andk-level of the level set function can extract the OD and OC contours accurately.


## 3. Numerical Calculation

We can minimize the energy functionalE alternately. For the parameters c1, c2, c3, bx, α, a, b, θ, x, and y fixed, we use the standard gradient descent method to minimize the energy functional E with respect to ϕ and obtain the corresponding gradient descent flow equation:(14)∂ϕ∂t=λ1δϕe1−λ2δϕ−δϕ−ke2−λ3δϕ−ke3+v1δϕ+v2δϕ−kdiv∇ϕ∇ϕ+μ∇2ϕ−αdiv∇ϕ∇ϕ−2γ1δϕHϕx−Hϕ1x−2γ2δϕHϕ1x−Hϕ1x,where ei is representative of eix=∫Kρx−yIx−byci2dy, which can be rewritten as ei=Kρ∗1I2−2ciIby∗Kρ+ci2by2∗Kρ. δx is the derivative of the Heaviside function. div⋅ is the divergence operator.For the parametersϕ, bx, α, a, b, θ, x, and y fixed, the energy functional E is minimized with respect to c1, c2, and c3, which are given as follows:(15)c1=∫bx∗KρIM1ϕydy∫bx2∗KρM1ϕydy,(16)c2=∫bx∗KρIM2ϕydy∫bx2∗KρM2ϕydy,(17)c3=∫bx∗KρIM3ϕydy∫bx2∗KρM3ϕydy.For the parameterϕ fixed, the distance regularization term Dϕ,α is minimized by using the gradient descent flow method. Therefore, the evolution equation with respect to α can be obtained:(18)∂αx∂t=μ1∇ϕx−αx+2μ2∇2αx,For the parametersϕ, c1, c2, c3, α, a, b, θ, x, and y fixed, we can compute the bias field bx by minimizing the energy functional E, which can expressed as(19)bx=∑i=13Kρ∗IMiϕxci∑i=13Kρ∗Miϕxci2.For the parametersϕ, c1, c2, c3, bx, and α fixed, the evolution function of a, b, θ, x, and y can be obtained by minimizing the energy functional E, which is given by(20)da1dt=2γ1∫ΩHϕ−Hϕ1δϕ1A12a13dx,(21)db1dt=2γ1∫ΩHϕ−Hϕ1δϕ1B12b13dx,(22)dx1dt=2γ1∫ΩHϕ−Hϕ1δϕ1A1cosθ1a12−B1sinθ1b12dx,(23)dy1dt=2γ1∫ΩHϕ−Hϕ1δϕ1A1sinθ1a12+B1cosθ1b12dx,(24)dθ1dt=2γ1∫ΩHϕ−Hϕ1δϕ1−A1B1a12+A1B1b12dx,(25)da2dt=2γ2∫ΩHϕ−Hϕ2δϕ2A22a23dx,(26)db2dt=2γ2∫ΩHϕ−Hϕ2δϕ2B22b23dx,(27)dx2dt=2γ2∫ΩHϕ−Hϕ2δϕ2A2cosθ2a22−B2sinθ2b22dx,(28)dy2dt=2γ2∫ΩHϕ−Hϕ2δϕ2A2sinθ2a22+B2cosθ2b22dx,(29)dθ2dt=2γ2∫ΩHϕ−Hϕ2δϕ2−A2B2a22+A2B2b22dx,where A1, A2, B1, and B2 are in the following forms:(30)A1=x−x1cosθ1+y−y1sinθ1,B1=−x−x1sinθ1+y−y1cosθ1,A2=x−x2cosθ2+y−y2sinθ2,B2=−x−x2sinθ2+y−y2cosθ2.In practical implementation, the regularized Heaviside functionHε is used to maintain the stability of numerical calculations, which is given by(31)Hε=121+2πarctanxε,and its derivative δεx is derived in the form below:(32)δεx=1πεε2+x2.In this paper, we use an explicit scheme to solve the numerical problem. All the spatial partial derivatives and∂ϕ/∂x and ∂ϕ/∂y are approximated by the central difference, and the temporal partial derivative ∂ϕ/∂t is approximated by the forward difference [34]. Then, equation (13) can be discretized as(33)ϕn+1=ϕn+Δtλ1δϕne1−λ2δϕn−δϕn−ke2−λ3δϕn−ke3+v1δϕn+v2δϕn−kdiv∇ϕn∇ϕn+μ∇2ϕn−αdiv∇ϕn∇ϕn−2γ1δϕnHϕnx−Hϕ1nx−2γ2δϕnHϕnx−Hϕ2nx,where Δt is the time step and n 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:(34)ϕn+1=ϕn+Δt1⋅Δϕn.The specific steps of the proposed method are presented as follows:Step 1: initialize level set functionsϕx, ϕ1x, and ϕ2x as signed distance functions; initialize the function αx and bias field bx.Step 2: calculatec1, c2, and c3 according to equations (15)–(17).Step 3: updatea, b, θ, x, and y according to equations (20)–(29).Step 4: evolve the level set functionϕx according to equation (14).Step 5: smooth the level set functionϕx using equation (34).Step 6: calculateαx using equation (18).Step 7: calculatebx 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ϕx converges or not, which is defined by(35)Q=LC0n−LC0n−1+LCkn−LCkn−1,whereQ is the sum of length changes of 0-level contour and k-level contour between n and n−1 iterations. When Q<τ (τ 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: Δt=0.1, Δt1=0.01, ρ=3, μ=1, μ1=1, μ2=0.5, λ1=1, λ2=3, λ3=1, γ1=100, γ2=100, v1=0.06∗255∗255, v2=0.01∗255∗255, and α=1. 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(36)Dice=2×AreaX∩YAreaX+AreaY.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 Figure3. 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.Figure 3
The comparison of the evolution process of level set functions in different methods. (a) The experimental results of the proposed method without prior elliptic shape constraint term. (b) The experimental results of the proposed method. From top to bottom in the figure are 100 iterations, 200 iterations, 300 iterations, and final results.
(a)
(b)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 Figures4–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.Figure 4
Experimental results of the OD and OC segmentation with mild glaucoma. (a) Initialization. (b) Multiphase C-V. (c) DR2LS. (d) DRLSEC. (e) Ground truth.
(a)
(b)
(c)
(d)
(e)Figure 5
Experimental results of the OD and OC segmentation with moderate glaucoma. (a) Initialization. (b) Multiphase C-V. (c) DR2LS. (d) DRLSEC. (e) Ground truth.
(a)
(b)
(c)
(d)
(e)Figure 6
Experimental results of the OD and OC segmentation with severe glaucoma. (a) Initialization. (b) Multiphase C-V. (c) DR2LS. (d) DRLSEC. (e) Ground truth.
(a)
(b)
(c)
(d)
(e)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 Table1. 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.Table 1
Comparison of segmentation performance with different methods.
Method
Average overlapping ratio of OD (%)
Average overlapping ratio of OC (%)
Multiphase C-V method
62.33
60.04
DR2LS
63.21
61.56
DR2LS with EC
68.82
63.48
DRLSEC without EC
67.70
64.96
DRLSEC
72.16
66.45


## 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.



---

*Source: 4836296-2019-10-24.xml*