Vessel bend-based cup segmentation in retinal images
Gopal Datt Joshi, Jayanthi Sivaswamy
CVIT, IIIT Hyderabad, India
gopal@research, jsivaswamy@iiit.ac.in
Kundun Karan, Prashanth R., S. R. Krishnadas
Aravind Eye Hospital, Madurai, India
krishnadas@aravind.org
Abstract
In this paper, we present a method for cup boundary
detection from monocular colour fundus image to help
quantify cup changes. The method is based on anatomi-
cal evidence such as vessel bends at cup boundary, con-
sidered relevant by glaucoma experts. Vessels are mod-
eled and detected in a curvature space to better handle
inter-image variations. Bends in a vessel are robustly
detected using a region of support concept, which au-
tomatically selects the right scale for analysis. A reli-
able subset called r-bends is derived using a multi-stage
strategy and a local spline fitting is used to obtain the
desired cup boundary. The method has been success-
fully tested on 133 images comprising 32 normal and
101 glaucomatous images against three glaucoma ex-
perts. The proposed method shows high sensitivity in
cup to disk ratio-based glaucoma detection and local
assessment of the detected cup boundary shows good
consensus with the expert markings.
1. Introduction
Early detection and treatment of retinal diseases are
crucial to avoid preventable vision loss. Glaucoma is
one of the most common causes of the blindness and
about 79 million in the world are likely to be afflicted
with glaucoma by the year 2020. It leads to irreversible
vision loss due to significant loss of optic nerve fibers.
Retinal nerve fibers converge to the optic disk (OD)
and form a cup-shaped region known as the cup. En-
largement of this cup with respect to OD is an important
indicator of glaucoma progression and hence ophthal-
mologists manually examine the OD and cup for eval-
uation. An automatic assessment of cup region from
colour fundus image (CFI) could reduce the workload
of experts and aid objective detection of glaucoma.
A CFI is a projection of retinal structures on 2D color
plane where OD appears as a bright circular or elliptic
region partially occluded by blood vessels as shown in
fig 1(a). Given the loss of depth information in a CFI,
glaucoma experts use a change in vessel morphology
Fig. 1: a) OD-centric CFI b) Cup boundary through r-bends.
as a reliable visual clue for determining cup boundary.
Thus cup segmentation is a time-consuming and chal-
lenging task. Figure 1(b) shows cup boundary marked
a glaucoma expert using vessel bends (highlighted by
arrows) information.
Much of existing work has mainly focused on OD
segmentation with only few attempts at cup segmenta-
tion. Cup boundary detection has been proposed using
3-D depth information [5]. Since 3-D images are not
easily available, Liu et. al. [3] proposed a cup bound-
ary estimation method using monocular CFI. In this ap-
proach, a potential set of pixels belonging to cup region
is first derived based on the reference colour obtained
from a manually selected point. Next, an ellipse is fit to
this set of pixels to estimate the cup boundary. A variant
of this method obtains the cup pixels via thresholding of
the green colour plane [3]. Cup boundary obtained via
ellipse fitting yields only coarse cup boundary. Further-
more, fixed thresholding is also not adequate to han-
dle large inter-image intensity variations that arise due
to complex imaging and physiological difference across
patients.
In order to address these problems, additional in-
formation such as small vessel bends (’kinks’) which
anatomically mark the cup boundary have been used in
[4]. Here, image patches are extracted around a esti-
mated cup boundary obtained in [3] and vessel pixels
are identified using edge and wavelet transform infor-
mation. Next, vessel bends, characterized by points of
direction change in the vessel pixels are found and used
Fig. 2: The proposed method
to obtain the cup boundary. This method is highly de-
pendent on the preliminary cup boundary obtained from
[3]. Furthermore, the statistical rule for selecting vessel
pixels is very sensitive to the inter-image variations.
Both appearance and anatomical knowledge are used
by the glaucoma expert to determine cup boundary in
different cup regions. Hence, we propose a method
that integrates both information under a common frame-
work. The cup is modeled as a region enclosing pallor
region (shown in fig. 1(a)) and defined by a boundary
passing through a sparse set of vessel bends. In the next
section, we will explain proposed method in detail.
2. Proposed Method
Our objective is to segment the cup region by us-
ing anatomical evidences considered relevant by the ex-
perts. As seen in fig. 1(b) and cyan points in fig. 3(b),
vessel bends can occur at many places within the OD re-
gion. However, only a subset of these points define the
cup boundary. We refer to this as relevant vessel bend or
r-bend. The first problem at hand is to find this subset.
We use multiple source of information for this purpose:
the pallor region which spatially defines the inner limit
of r-bend, bending angle and location in the OD region.
A second problem is that the anatomy of the OD region
is such that the r-bends are non-uniformly distributed
across a cup boundary with more points on the top and
bottom; they are mostly absent in the nasal side and
very few in number in the temporal side. We propose
a local interpolating spline to naturally approximate the
cup boundary in regions where r-bends are absent. Fig-
ure 2 shows an overview of the proposed method. The
energy minimization based deformable models are not
adequate due to absence of certain edge or region based
information associated with the cup region to derive an
energy functional.
Optic Disk Segmentation: The OD region is lo-
calised from red channel of the CFI and morphological-
based pre-processing step is performed to reduce the ef-
fect of vessels within OD. We use the region-based ac-
tive contour approach in [1] for segmenting the OD. The
contour C(s) : [0, 1] → IR
2
, which is a piecewise pa-
rameterized C
1
curve, is evolved using an energy func-
tional defined as
F (c
+
, c
−
, C) = µ . Length(C) (1)
+λ
+
Z
inside(C)
|I
0
(x, y) − c
+
|
2
dxdy
+λ
−
Z
outside(C)
|I
0
(x, y) − c
−
|
2
dxdy
Fig. 3: a) Angle of a vessel bend, b) uniform pallor sam-
ples(red), bend points(cyan) and c) fitted circle(red) and po-
tential r-bends.
where c
+
and c
−
are unknown constants representing
the average value of I
0
inside and outside the curve, re-
spectively. The parameters µ ≥ 0 and λ
+
, λ
−
≥ 0;
are weights for the regularizing and the fitting terms,
respectively. For curve evolution, the level set formula-
tion is used where the motion is governed by mean cur-
vature [1]. For further processing steps, we only con-
sider green colour plane of segmented OD region.
Medial axis detection: The OD region has both thick
and thin vessels and due to large inter-image variation,
detecting both kinds is difficult. Hence, we embed the
image in a 3D space and formulate the blood vessel de-
tection as a problem of trench detection in the intensity
surface. The selection of this space gives robustness to
the image variations and detection is solely driven by
trench shape and directional continuity associated with
a vessel structure. Trenches are regions characterized
by high curvature, oriented in a particular direction.
The curvature is computed using surface tangent
derivative [2] defined as: Υ(x) =
d
2
y/dx
2
1+(dy/dx)
2
. For each
point, Υ is computed in 4 different directions. The max-
imum value of the responses Υ
max
and corresponding
orientation α (perpendicular to the vessel direction) are
retained and further assessed to obtain trench points. A
point is declared as a trench if value of Υ
max
is greater
than both threshold value t and the values of neighbor-
ing pixels in α direction.
For the robust detection of low contrast vessels, we
employ a two-phase thresholding scheme in which first,
a high value of t is applied to get high contrast vessel
points (set-1). Then, low value of t is applied to get
a new set of low contrast vessel points (set-2). Points
in set-2 which are found connected to the set-1 are in-
cluded in the final set along with set-1. This strat-
egy helps in successfully extracting low contrast vessels
while rejecting noise. The final trench points give a me-
dial axis based representation of vessel structure which
is more precise in quantifying vessel bends compared
to edge-based representation. The next task is to extract
vessel bends from this representation.
Vessel Bend detection: The amount of bending in
vessels varies according to the caliber of vessel. Thin
vessels show significant bending compared to a thick
vessel. This is due to the fact that thick vessels are
more rigid. The selection of appropriate scale for de-
tecting both types of bend is crucial because bend in a
thick vessel is apparent only at a larger scale compared
to a bend in thin vessel. We employ a dynamic region
of support (ROS) based scheme to find the appropriate
scale to analyse a candidate point.
First, we extract vessel segments terminated by end
and/or junction points. For each segment, we compute
1D shape (curvature) profile and locate the local max-
ima. These local maxima constitute a candidate set of
bends b. A ROS for any b
i
is defined as a segment of
vessel around b
i
and bound on either side by the nearest
curvature minimum. Choosing the bounds to be based
on curvature minima automatically ensures the size of
the ROS to be large for thick vessels and small for thin
vessels. The angle of bend θ is then computed as the
angle between the lines joining a bend point and the
centers of mass on both sides of the ROS. The center of
mass of an arm is defined by the mean position of pix-
els on the arm (illustrated in fig. 3(a)). Since only ves-
sels bending into the cup are of interest, bends above
θ = 170
◦
are eliminated from the candidate set. The
detected vessel bends in a sample image are highlighted
in fig. 3(b) with cyan markers.
Multi-stage selection of r-bends The task of identi-
fying the r-bends from b
i
is performed in two stages,
to reduce the required analysis, by utilizing anatomical
knowledge associated with r-bends. In the first stage, a
coarse selection is done based on a bend’s proximity to
the pallor region. In the second stage, the spatial posi-
tion and bending information are used to identify the set
of r-bends.
First stage: Let p : (x
p
, y
p
) be a set of points within
the pallor region. These were found by retaining the
top 25% of the bright pixels within the OD. Next, let
b : (x
b
, y
b
) be the locations of the bends b
i
. The re-
gion containing potential r-bends is localised by find-
ing a best-fit circle (in least-square sense) to the set
of points (x, y) = {p, b}. Let the circle have center
(x
c
, y
c
) and radius R. The minimisation of error func-
tion S =
P
i
((x
i
− ¯x) − x
c
)
2
+ ((y
i
− ¯y) − y
c
)
2
− R
2
gives the unknown parameters x
c
, y
c
and R, where,
¯x =
1
n
P
i
x
i
, ¯y =
1
n
P
i
y
i
and n number of points.
The bends which lie in the vicinity of this circle (inside
and outside) are passed to the next stage. Figure 3(c)
shows sample candidate r-bends obtained in this stage.
Second stage: Each candidate bend is analysed in
terms of its sector-wise location (as in fig. 3(c)) and
its parent vessel orientation. This analysis is based on
anatomical knowledge that bends formed by vertical
vessels in sec-1&3 and horizontal vessels in sec-2&4
are the probable r-bends. The final refined set of r-bends
Fig. 4: a) Estimated cup boundary, b) final OD and cup bound-
ary.
Fig. 5: Detected cup boundary.
is found as follows: A sector is radially analysed with
a step size of 20
◦
and in each step, only bends formed
by vessels with the
0
correct
0
orientation are retained. If
multiple bends remain, then the bend with smaller value
of theta θ is selected as thin, rather than thick, vessel
bends are more reliable indicators for the cup boundary.
These usually occur in the diagonal region between two
sectors.
2D spline interpolation Typically, r-bends are sparse
and not uniformly distributed across the sectors. In their
absence, experts use their clinical knowledge (experi-
ence of direct 3D cup examination) to approximate a
cup boundary. Hence, it is difficult to get the cup bound-
ary in the regions with no r-bends. We choose a lo-
cal cubic cardinal spline, which is a generalisation of
Catmull-Rom spline, with a shape parameter t. The pa-
rameter t helps control the bending behaviour and thus
the shape according to the sector. The value of t is kept
high in sectors 2&4 as they usually have low vessel den-
sity (r-bends) compared to sector 1&3. A closed-form
2D spline curve is obtained by considering, sequen-
tially, a subset of r-bends. Figure 4(a) shows the inter-
polated cup boundary passing through the r-bends and
Fig. 4(b) shows final obtained boundaries for a sample
OD region.
3. Experimentation Results
The proposed method was evaluated on a dataset of
retinal images collected from an ongoing pilot study in
collaboration with a local eye hospital. The dataset has
32 normal and 101 glaucomatous (total of 133) images.
All images were taken under a fixed protocol with 30-
degree field of view, centered on the OD. For each im-
age, ground truth was collected from three glaucoma
experts, referred to as E1, E2 and E3 with experience of
3, 5 and 20 years, respectively.
Figure 5 shows the detected cup boundary against
Expert-1 Expert-2 Expert-3
Cat/No. µ σ µ σ µ σ
N/32 0.56 0.79 0.23 0.19 0.18 0.28
G/101 0.06 0.20 0.05 0.35 0.03 0.34
Total/133 0.09 0.80 0.09 0.33 0.01 0.34
Table 1: Mean µ and σ in CDR estimation.
three experts on a sample image (overlaid on segmented
OD region). The proposed method successfully detects
r-bends formed on both thick and thin vessel. The cup
boundary at r-bends is closer to the experts’ marked
boundaries (fig. 5(a-2)), whereas in regions where they
are absent the interpolated result is unable to match the
boundary marked by the experts. We observed some
challenging situations where our detected r-bends are
not considered relevant by experts. For instance, in fig.
5(a-1) boundaries marked by experts are away from the
detected r-bends though there was not 2D clue present
to support their markings. These suggest the role of
prior 3D knowledge of the cup being used by an expert
to determine cup boundary.
Traditionally, the cup-to-disk ratio(CDR) in the ver-
tical direction is used to quantify cup enlargement. This
was computed from the obtained boundaries and com-
pared against that from each of the experts. The average
mean µ and standard deviation σ of the CDR error for
32 normal and 101 glaucomatous images is shown in
Table. 1. The method gives less estimation error against
expert-3. The average error µ, σ over all three experts
is 0.323, 0.420 for normal and 0.046, 0.296 glaucoma-
tous images. These figures show that high accuracy
of method in estimating CDR for glaucomatous images
compared to normal image indicating high sensitivity in
glaucoma detection.
The CDR measure is inadequate to assess local cup
changes which is of clinical interest. We examined this
by computing sector-wise cup segmentation accuracy.
An area (pixel) overlap-based method was used to
compute the precision(P)-recall(R) measures to assess
the cup segmentation in each of the 4 sectors. Table 2
shows the sector-wise P and R measure computed over
133 images. These figures indicate that the method
is consistent across four sectors and with an average
precision(0.79) and recall(0.87).
In order to gain some insight into the reasons behind
this performance we analysed the inter observer vari-
ance among the experts by taking E3 with 20 years ex-
perience, as our gold standard. The difference in cup
radius between the gold standard and the experts were
found for each point on the boundary (defined by the
angle α) for all 133 images. These are plotted in fig.
6. It can be seen that on average the variance among
the experts is E1 = 24.06 and E2 = 23.91. This indi-
cates a fair degree of disagreement between experts and
Sec-1 Sec-2 Sec-3 Sec-4
P R P R P R P R
E-1 0.81 0.81 0.74 0.87 0.77 0.89 0.82 0.84
E-2 0.87 0.81 0.82 0.83 0.82 0.85 0.85 0.85
E-3 0.75 0.89 0.72 0.89 0.69 0.94 0.72 0.92
Avg 0.81 0.84 0.76 0.86 0.76 0.89 0.80 0.87
Table 2: Sector-wise precision and recall measures for cup
segmentation.
Fig. 6: Angular cup radius assessment against expert-3.
attests to the complexity of the problem.
4 Discussion and Conclusion
In this paper, we presented a novel cup boundary
detection method using r-bends information. Trench
based vessel modeling and ROS-based bend detection
that have been employed result in robustness to vary-
ing thickness of the vessels. Final cup boundary is ob-
tained using applying local spline interpolation on the
detected r-bends. Assessment results show that our
method matches quite well with the most experienced
expert’s assessment. At a local level, the segmentation
P/R figures are 0.79/0.87. It is observed that in the re-
gions with no certain 2D clues, there is less consensus
on the cup boundary between our method and experts
and also within experts. This signals the ambiguity in
2D information and the importance of 3D information
in cup segmentation which will be investigated in our
future work.
References
[1] T. Chan and L. Vese. Active contours without edges.
IEEE Trans. Image Processing, 10(2):266–277, 2001.
[2] S. Garg, J. Sivaswamy, and S. Chandra. Unsupervised
curvature-based retinal vessel segmentation. Proc. ISBI,
pages 344–347, 2007.
[3] J. Liu, D. Wong, J. Lim, H. Li, N. Tan, and T. Wong.
Argali- an automatic cup-to-disc ratio measurement sys-
tem for glaucoma detection and analysis framework.
Proc. SPIE, Medical Imaging, pages 72603K–8, 2009.
[4] D. Wong, J. Liu, J. H. Lim, H. Li, X. Jia, F. Yin, and
T. Wong. Automated detection of kinks from blood ves-
sels for optic cup segmentation in retinal images. Proc.
SPIE, Medical Imaging, page 72601J, 2009.
[5] J. Xu, O. Chutatape, E. Sung, C. Zheng, and P. Chew. Op-
tic disk feature extraction via modified deformable model
technique for glaucoma analysis. Pattern Recognition,
40(7):2063–2076, 2007.
