Level Set Segmentation of Retinal OCT Images

Bashir Isa Dodo

, Yongmin Li

, XiaoHui Liu

and Muhammad Isa Dodo

Department of Computer Science, Brunel University, U.K.

Katsina State Institute of Technology and Management, Katsina, Nigeria

Keywords:

Image Segmentation, Level Set, Evolution Constrained Optimisation, Optical Coherence Tomography.

Abstract:

Optical coherence tomography (OCT) yields high-resolution images of the retina. Reliable identiﬁcation of

the retinal layers is necessary for the extraction of clinically useful information used for tracking the progress

of medication and diagnosis of various ocular diseases. Many automatic methods have been proposed to aid

with the analysis of retinal layers, mainly, due to the complexity of retinal structures, the cumbersomeness of

manual segmentation and variation from one specialist to the other. However, a common drawback suffered by

existing methods is the challenge of dealing with image artefacts and inhomogeneity in pathological structures.

In this paper, we embed prior knowledge of the retinal architecture derived from the gradient information, into

the level set method to segment seven (7) layers of the retina. Mainly, we start by establishing the region

of interest (ROI).The gradient edges obtained from the ROI are used to initialise curves for the layers, and

the layer topology is used in constraining the evolution process towards the actual layer boundaries based on

image forces. Experimental results show our method obtains curves that are close to the manual layers labelled

by experts.

1 INTRODUCTION

Optical coherence tomography ﬁrst introduced by

(Huang et al., 1991) is a noninvasive imaging tech-

nique that provides cross-sectional images of the

retina with an acquisition speed of approximately

25,000 A-scans per second, and an axial image reso-

lution of approximately 5-7µm (Raftopoulos and Trip,

2012; Jaffe, 2012; Adhi and Duker, 2013). In cur-

rent ophthalmology, identifying various layers of the

retina on OCT has become a vital tool for diagnos-

ing and tracking the progress of medication of various

visual impairments (Adhi and Duker, 2013). Since

the manual segmentation is not only tedious but also

subjected to intra- and inter-grader variance (Yazdan-

panah et al., 2011), many automatic methods have

been proposed to assist with the segmentation pro-

cess(Sun et al., 2016). Current research in retinal

OCT segmentation focuses on improving various as-

pects ranging from computational time, the number of

layers segmented, use of prior knowledge and compu-

tational complexity to mention a few. In general, seg-

mentation is partitioning based on some image char-

acteristics. How these characteristics are deﬁned de-

termines the computation burden of the algorithm.

In some cases, this computational burden is reduced

through dynamic programming (Chiu et al., 2010) or

topology modiﬁcation (Liu et al., 2018). The number

of questions (conditions) an algorithm has to check or

satisfy is usually the most crucial factor.

Particularly, Markov Boundary Model

(Koozekanani et al., 2001), later extended by

(Boyer et al., 2006), geodesic distance (Duan et al.,

2017), level sets (Wang et al., 2015; Novosel et al.,

2013), graph-based methods (Chiu et al., 2010;

Garvin, 2008; Dodo et al., 2017; Dodo et al., 2018),

and machine learning (Lang et al., 2013) have been

used in the segmentation of retinal OCT images.

Although the level sets method has automatic topo-

logical handling, the steps can be computationally

expensive (Shi and Karl, 2005), while adding

complex constraints in the segmentation method

usually increases the complexity of an algorithm.

In this work, we incorporate simple yet efﬁcient

topological constraints to the evolution process

of the level set method, to improve accuracy and

reduce computational complexity for OCT image

segmentation.

The method we propose is based on the fol-

lowing considerations: 1. The image gradients are

used to initialise curves in order to handle under-

segmentation and over-segmentation of the image; 2.

Dodo, B., Li, Y., Liu, X. and Dodo, M.

Level Set Segmentation of Retinal OCT Images.

DOI: 10.5220/0007577600490056

In Proceedings of the 12th International Joint Conference on Biomedical Engineering Systems and Technologies (BIOSTEC 2019), pages 49-56

ISBN: 978-989-758-353-7

The evolution of a curve is explicitly based on layer

arrangements and implicitly based on OCT topology.

This means for each image the initial contours are

speciﬁc to the image under investigation, while the

forces in the normal direction and the topology con-

strains guide the contours towards layer boundaries.

Our method segments an OCT image into 7 seg-

ments, relating to: Nerve Fibre Layer (NFL); Gan-

glion Cell Layer + Inner Plexiform Layer + In-

ner Nuclear Layer (GCL+IPL+INL); Outer Plexiform

Layer (OPL); Outer Nuclear Layer to Inner Segment

(ONL+IS), Outer Segment (OS) and Retinal Pigment

Epithelium (RPE). Locations of these layers on an

OCT image are shown in Figure 1. The rest of the

paper is organised as follows. Section 2 discusses the

proposed method in details. Experimental results and

discussions are treated in Section 3. Finally, conclu-

sions are drawn in Section 4.

2 PROPOSED METHOD

This section describes our approach of segmenting

retinal OCT images. A schematic representation of

the method is illustrated in Figure 2, and details of

each step are elaborated in the subsequent subsec-

tions.

2.1 Pre-processing

The pre-processing steps are illustrated in Figure 3.

Each OCT B-scan image I is ﬁrst enhanced with a

Gaussian ﬁlter to reduce the image noise. The lay-

ers targeted in our study lies within the total retinal

thickness (TRT), which starts from the Internal Limit-

ing Membrane(ILM) to the posterior boundary of the

Retinal pigment epithelium (RPE), i.e. the boundary

between the retinal nerve ﬁbre layer and the vitreous,

and the boundary between the RPE and the choroid

regions respectively. It is commonly accepted that the

NFL, IS-OS and RPE exhibits high reﬂectivity in an

OCT image (Chiu et al., 2010; Lu et al., 2011; Tian

et al., 2015), based on experiments the ILM and RPE

exibits the highest transitions from dark-bright and

bright-dark, respectively (Dodo et al., 2018). Based

on this understanding of the retinal structure, the ILM

and RPE are identiﬁed using the shortest path (Dijk-

stra, 1959), by searching for the highest transitions on

two separate adjacency matrices(Chiu et al., 2010).

Using the ILM and RPE points the image is

cropped to I

cropped

, such that only the Region of in-

terest (ROI) with useful layer information is remain-

ing. This helps in dealing with layer-like structures

outside the ROI and the computational cost associ-

ated with handling image background in segmenta-

tion. The pre-processing is vital in our segmentation

process because only the actual layer properties im-

pact the evolution of the curve. The next operation on

the image is to generate a mask I

mask

of the cropped

image, and then multiply it by the original image I.

The examples of resulting images from this step are

shown in ﬁgure 3, column 4, and expressed by the

equation below:

processed

= I

mask

∗ I (1)

The process in this subsection is essential because

only the layer structures are obtained when the gradi-

ent of the image is acquired. One of the signiﬁcant

roles of the pre-processing is to eliminate the need for

handling background as depicted in ﬁgure 4 and fur-

ther discussed in the next few subsections 2.2 and 2.3.

From Figure 4 (a) the background and image noise

will affect the segmentation processes (because the

area highlighted in red will be initialised), and for this

reason we establish the ROI only based on part of the

image that contains useful information Figure 4 (b).

The estabilishment of the ROI also complements the

thresholding and reﬁnement processes in the layer ini-

tialisation stage. We use the size of the cropped image

to reduce computational time further and to eliminate

the need for storing idle points.

2.2 Boundary Initialisation

To initialise contours we obtain the vertical gradient

∇I

processed

of the processed image and threshold it by

a constant T , in our case T = 0.0018. The value of

T should ideally be low, to avoid getting more com-

ponents in the GCL - IPL regions which will nega-

tively impact the segmentation results. We obtain the

edges of the thresholded gradient T G image 5(a), and

then reﬁne it in two simple steps: ﬁrst, by area open-

ing ((Vincent, 1994), where any component less than

P pixels (P = 30pixels) is deleted to remove small

objects from the image (most especially the GCL to

IPL regions); second, we extrapolate incomplete layer

lines to span the image horizontally. This is carried

out by linking broken lines to the closet neighbouring

points, where each line starts from the ﬁrst column

and ends at the last column of the image, such that

only complete layers are initialised, ﬁgure 5(b). The

initialisation and evolution of complete layer lines ex-

clusively is an essential factor, which further ensures

accurate segmentation, i.e. no merging or splitting of

boundaries is allowed, which prevents under- or over-

segmentation. Without losing context, the reﬁnement

step can be ignored, if the layers from the GCL to INL

BIOIMAGING 2019 - 6th International Conference on Bioimaging

Figure 1: Location of Nerve Fibre Layer (NFL); Ganglion Cell Layer + Inner Plexiform Layer + Inner Nuclear Layer

(GCL+IPL+INL); Outer Plexiform Layer (OPL); Outer Nuclear Layer to Inner Segment (ONL+IS), Outer Segment (OS),

Retinal Pigment Epithelium (RPE) and the total retinal thickness, on an OCT image.

Figure 2: Schematic representation of the proposed level set approach.

are targets of the method. However, this will require

a condition for handling the splitting and merging of

boundaries or an alternate measure to correctly iden-

tify which layers are segmented. Moving further, the

edges of the reﬁned image serve as initial curves such

that the number of identiﬁed regions in the ﬁnal out-

put cannot exceed the number of the initial curves.

An important point to note is that, if T is reduced,

then P is to be increased, mainly because the size of

the small objects in 5 (b) will increase with a smaller

Level Set Segmentation of Retinal OCT Images

Figure 3: Preprocessing steps showing: Column 1 - Enhanced images; Column 2 - identiﬁed ILM (red) and RPE (Green);

Column 3 - image masks I

mask

; and Column 4 - Cropped images I

cropped

. Row 1 - Nasal region; Row 2 - Foveal Region; and

Row 3- Temporal Region.

Figure 4: Gradient of full image ( Figure - 3 Column 1, Row 1) with background noise and layer-like structures in red (a), and

thresholded gradient of preprocessed (Figure 3 - Column 4, Row 1) image TG with ROI only (b).

T . Each boundary curve C

is therefore represented

by a collection of C

(x,y) on the image.

2.3 Segmentation of OCT

Next, each initial boundary curve C

is evolved de-

pending on a speed ﬁeld F based on the following

differential equation (Shi and Karl, 2005) :

= F

N (2)

Where

N is the normal of the curve pointing out-

ward. The speed ﬁeld F is made of an external speed

derived from the image data and a characteristic speed

BIOIMAGING 2019 - 6th International Conference on Bioimaging

Figure 5: Edges before reﬁnement (a), and reﬁned edges used for contour initialisation (b).

based on C

. We associate F and the ensuing evolu-

tion with a gradient descent solution to solve the min-

imisation problem based on a Mumford-shah model

evolution perspective (Tsai et al., 2001). This means

the curve C

will evolve until it gets to a local minima

bmin

of the energy, i.e. static point of the dynamic

equation (2). Adapting from (Shi and Karl, 2005), we

represent a layer boundary uniquely through two lists

of inside L

and outside L

out

points of C

. Which are

deﬁned as:

out

= {x|φ(x) > 0 and ∃y ∈ N(x) : φ(y) < 0}

= {x|φ(x) < 0 and ∃y ∈ N(x): φ(y) > 0}

where N(x) is a distinct neighbourhood of x, in the

level set function φ at pixel x. Based on this deﬁnition,

a positive force switches a point from L

to L

out

and a

negative force switches a point from L

out

to L

. Each

point (x, y) in level set function is deﬁned in relation

to the curve C

as follows (Shi and Karl, 2005):

φ =











3, if x,y is outside C

and x,y 6∈ L

out

;

1, if x,y ∈ L

out

;

−3, if x,y is inside C

and x,y 6∈ L

;

−1, if x,y ∈ L

(3)

Based on the deﬁnition of φ in equation (3), it

is undemanding to recognise the location of a point

(x,y) on the image in relation to C

. In our formu-

lation, we use a 2D list to represent initial bound-

ary points, and to save the positions of ﬁnal bound-

ary points for straightforward mapping in generating

the ﬁnal image output. Alternatively, a 1-D list can

be used to save the boundaries point at position φ, as

suggested in (Liu et al., 2018). A boundary position

(x,y) of φ, can either expand or shrink based on:











Expand(x,y) : C

(x,y) := C

(x,y) + 1

Shrink(x,y) : C

(x,y) := C

(x,y) − 1

(4)

The evolution of each point is inﬂuenced by the

image forces computed by a fast gradient vector

ﬁeld(Wang and Boyer, 2012) and the topology con-

straints to be described in the next subsection 2.4.

2.4 Topology Constrains

As highlighted earlier in subsection 2.1 the ordering

of the layers must be preserved. Taking into account

the architecture of the OCT image, boundary C

always below C

for any given boundary points C

and C

, i.e. a point (x,y) on the curve will neither

Shrink nor Expand if it makes C

(x,y) ≤ C

(x,y).

Hence, with our appropriate initialisation, we enforce

the topology requirement by carrying out this simple

topology validation before either shrinking or expand-

ing a boundary. Finally, we employ an intuitive ap-

proach to ensure this topology is preserved, by addi-

tionally reﬁning the topology constraint in the vertical

direction:

1. Because each layer boundary spans the image hor-

izontally (one boundary point per column) we add

a condition for evolving a boundary point C

(x,y)

to a new boundary point C

bmin

(x,y). We restrict

Expand(x, y) if its neighbour points are u con-

secutive points above it; do not Shrink(x, y) if its

neighbour points are u consecutive points below

it;

Level Set Segmentation of Retinal OCT Images

2. Looking at the sample of initial layer boundaries

in ﬁgure 5, a boundary point C

(x,y) is limited

to a maximum of v operations (either Expand or

Shrink) consecutively in the vertical direction).

The parameters u and v are two prior constants, in

our experiments u and v are set to 3 and 20 respec-

tively. The parameter u aids with boundary smooth-

ness and avoiding peaks for Expand(x, y) or valleys

for Shrink(x, y) on the boundaries, while v further en-

sures the layered architecture is preserved. Addition-

ally, this is why our layer initialisation is ideal be-

cause the starting points are based on the individual

image. The topology constraints facilitate the evolu-

tion because the validation is performed before ex-

panding or shrinking a boundary. Perhaps, this might

not be ideal for abnormal structures. However, con-

sidering the ordering of the layers where C

2 will al-

ways be below C

1 the layers will move together even

in the case of abnormal retinal structure. The pseudo-

code of our algorithm is illustrated in Algorithm 1.

3 RESULTS AND DISCUSSIONS

We applied our method to 200 macula OCT images.

The original size of each image is 512 by 992 pixels,

with a resolution of 16 bits per pixel. We crop the

image in the pre-processing stage to improve the re-

sults of our method. Comparison to the ground truth

labelling is carried out on full image size. In our ex-

periments N(x) = 8 neighbourhood, mainly, because

only the layers are remaining in the cropped image

and the effect of inhomogeneity is reduced. Experi-

mental results show that our method successfully seg-

ments seven (7) layers of the retina. Samples of the

method output are shown in ﬁgure 6.

Table 1 shows the mean and standard deviation

for the performance of the method compared to the

labelling of manual graders. The values show the

promising performance of our method in converging

at curves C

bmin

very close to the actual layer bound-

aries. Notably, the RNL thickness is used for diag-

nosing major eye diseases such as glaucoma, and the

mean (0.951) and standard deviation (±0.022) of dice

coefﬁcient for this layer is reassuring.

Moreover, it can be deduced that our method is

consistent in identifying the layer boundaries from

the distribution of the values in ﬁgure 7. Consider-

ing the second quarterlies of the NFL, IS and RPE

begin at ≥ 0.900 further attests to the optimum per-

formance of our method, except for few instances in

the GCL+IPL+INL and OPL layers, where the dice

score is below 0.800. However, in few instances the

method could not properly identify the GCL-INL and

Figure 6: From top to Bottom: Sample results from Nasal,

Foveal and Temporal regions respectively.

Table 1: Performance evaluation: Mean and Standard Devi-

ation (STDEV) of Dice Coefﬁcient on 200 B-Scan images

(Units in pixels).

RetinalLayer Mean ST DEV

NFL 0.951 ± 0.022

GCL+IPL+INL 0.879 ± 0.031

OPL 0.892 ± 0.032

ONL 0.907 ± 0.030

IS 0.932 ± 0.017

OS 0.920 ± 0.028

RPE 0.934 ± 0.021

the OPL due to some of the small components not

been removed. The method avoids over and under

segmentation, due to our layer initialisation and topo-

logical constraint, which prevents merging or splitting

of boundaries. On the other hand, the method will

come in short when compared to studies targeting the

choroidal region, which is believed to provide details

on some of the visual impairments (Sun et al., 2016).

BIOIMAGING 2019 - 6th International Conference on Bioimaging

Algorithm 1: Boundary Evolution.

1: Initialise Boundaries

2: loop:

3: if evolution will not make C

(x,y) ≤ C

(x,y) then

4: %% Shrink Boundary

5: if neighbours of C

(x,y) 6= consecutive u points below it then

6: if C

(x,y) has not moved v consecutive points in the vertical direction then

7: if force at point is negative then

8: Shrink(x, y)

9: %% Expand Boundaries.

10: if neighbours of C

(x,y) 6= consecutive u points above it then

11: if C

(x,y) has not moved v consecutive points in the vertical direction then

12: if force at point is positive then

13: Expand(x, y)

14: %% C

(x,y) are at relative local minima.

15: if no changes made to all C

(x,y) then

16: break

Figure 7: Box plot of mean Dice Coefﬁcient distribution for the seven (7) layers.

4 CONCLUSIONS

We have presented an automatic level set method

for retinal OCT segmentation. The proposed work

separates retinal OCT images into seven (7) non-

overlapping layers. Our approach has explored image

segmentation using level set from the point of initial-

isation, and the constraining of curve evolution based

on retinal layer topology explicitly. Reﬁned edges

of gradients images are used to initialise the curves.

Image forces constrained by the topological architec-

ture of the OCT are used to guide the evolution of

the curve to its minimum layer boundaries. These

two components ensure the boundaries obtained by

the method are close to the actual features of interest.

The proposed method takes advantage of handling the

obstruction of image background noise in the pre-

processing stage, consequently making the segmen-

tation process to suffer less from the image artefacts.

Additionally, reﬁnement of the gradient edge infor-

mation ensures only the targeted layers are initialised

at the beginning of the evolution process. Experimen-

tal results show that the proposed approach success-

fully segmented the target layers from OCT images,

and the segmentation results are close to the manu-

ally labelled ground-truth.Future work will seek to in-

clude the GCL to IPL and choroid regions in the ROI,

which has to do with implicitly assigning the param-

eters T and P based on the image as opposed to con-

stants used in our approach.

REFERENCES

Adhi, M. and Duker, J. S. (2013). Optical coherence

tomography–current and future applications. Current

Level Set Segmentation of Retinal OCT Images

opinion in ophthalmology, 24(3):213.

Boyer, K. L., Herzog, A., and Roberts, C. (2006). Auto-

matic recovery of the optic nervehead geometry in op-

tical coherence tomography. IEEE Transactions on

Medical Imaging, 25(5):553–570.

Chiu, S. J., Li, X. T., Nicholas, P., Toth, C. a., Izatt, J. a., and

Farsiu, S. (2010). Automatic segmentation of seven

retinal layers in SDOCT images congruent with expert

manual segmentation. Optics express, 18(18):19413–

19428.

Dijkstra, E. W. (1959). A note on two problems in connex-

ion with graphs. Numerische mathematik, 1(1):269–

271.

Dodo, B. I., Li, Y., Eltayef, K., and Liu, X. (2018). Graph-

cut segmentation of retinal layers from oct images. In

Proceedings of the 11th International Joint Confer-

ence on Biomedical Engineering Systems and Tech-

nologies - Volume 2: BIOIMAGING,, pages 35–42.

INSTICC, SciTePress.

Dodo, B. I., Li, Y., and Liu, X. (2017). Retinal oct im-

age segmentation using fuzzy histogram hyperboliza-

tion and continuous max-ﬂow. In 2017 IEEE 30th In-

ternational Symposium on Computer-Based Medical

Systems (CBMS), pages 745–750. IEEE.

Duan, J., Tench, C., Gottlob, I., Proudlock, F., and Bai,

L. (2017). Automated segmentation of retinal lay-

ers from optical coherence tomography images using

geodesic distance. Pattern Recognition, 72:158 – 175.

Garvin, M. K. (2008). Automated 3-D segmentation and

analysis of retinal optical coherence tomography im-

ages. PhD thesis - The University of Iowa.

Huang, D., Swanson, E. A., Lin, C. P., Schuman, J. S.,

Stinson, W. G., Chang, W., Hee, M. R., Flotte, T.,

Gregory, K., and Puliaﬁto, C. A. (1991). Optical

coherence tomography. Science (New York, N.Y.),

254(5035):1178–81.

Jaffe, G. J. (2012). OCT of the Macula: An expert provides

a primer on useful scans, identifying artifacts and time

domain vs. spectral domain technology. Reinal Physi-

cian, pages 10–12.

Koozekanani, D., Boyer, K., and Roberts, C. (2001). Reti-

nal thickness measurements from optical coherence

tomography using a Markov boundary model. IEEE

Transactions on Medical Imaging, 20(9):900–916.

Lang, A., Carass, A., Hauser, M., Sotirchos, E. S., Cal-

abresi, P. a., Ying, H. S., and Prince, J. L. (2013).

Retinal layer segmentation of macular OCT images

using boundary classiﬁcation. Biomedical optics ex-

press, 4(7):1133–52.

Liu, Y., Carass, A., Solomon, S. D., Saidha, S., Calabresi,

P. A., and Prince, J. L. (2018). Multi-layer fast level

set segmentation for macular oct. In 2018 IEEE

15th International Symposium on Biomedical Imaging

(ISBI 2018), pages 1445–1448.

Lu, S., Yim-liu, C., Lim, J. H., Leung, C. K.-s., and Wong,

T. Y. (2011). Automated layer segmentation of op-

tical coherence tomography images. Proceedings -

2011 4th International Conference on Biomedical En-

gineering and Informatics, BMEI 2011, 1(10):142–

146.

Novosel, J., Vermeer, K. A., Thepass, G., Lemij, H. G., and

Vliet, L. J. V. (2013). Loosely Coupled Level Sets

For Retinal Layer Segmentation In Optical Coherence

Tomography. IEEE 10th International Symposium on

Biomedical Imaging, pages 998–1001.

Raftopoulos, R. and Trip, A. (2012). The Application of

Optical Coherence Tomography ( OCT ) in Neurolog-

ical Disease. Advances In Clinincal Neuroscience and

Rehabilitation, 12(2):30–33.

Shi, Y. and Karl, W. C. (2005). A fast level set method

without solving pdes [image segmentation applica-

tions]. In Proceedings. (ICASSP ’05). IEEE Interna-

tional Conference on Acoustics, Speech, and Signal

Processing, 2005., volume 2, pages ii/97–ii100 Vol. 2.

Sun, Y., Zhang, T., Zhao, Y., and He, Y. (2016). 3d auto-

matic segmentation method for retinal optical coher-

ence tomography volume data using boundary surface

enhancement. Journal of Innovative Optical Health

Sciences, 9(02):1650008.

Tian, J., Varga, B., Somfai, G. M., Lee, W. H., Smiddy,

W. E., and DeBuc, D. C. (2015). Real-time automatic

segmentation of optical coherence tomography vol-

ume data of the macular region. PLoS ONE, 10(8):1–

20.

Tsai, A., Yezzi, A., and Willsky, A. S. (2001). Curve evolu-

tion implementation of the mumford-shah functional

for image segmentation, denoising, interpolation, and

magniﬁcation. IEEE Transactions on Image Process-

ing, 10(8):1169–1186.

Vincent, L. (1994). Morphological area openings and clos-

ings for grey-scale images. In Shape in Picture, pages

197–208. Springer.

Wang, C., Wang, Y., Kaba, D., Wang, Z., Liu, X., and Li., Y.

(2015). Automated layer segmentation of 3d macular

images using hybrid methods. In Proc. International

Conference on Image and Graphics. Tianjing, China.,

volume 9217, pages 614–628.

Wang, Q. and Boyer, K. L. (2012). The active geometric

shape model: A new robust deformable shape model

and its applications. Computer Vision and Image Un-

derstanding, 116(12):1178 – 1194.

Yazdanpanah, A., Hamarneh, G., Smith, B. R., and Sarunic,

M. V. (2011). Segmentation of intra-retinal layers

from optical coherence tomography images using an

active contour approach. IEEE Transactions on Med-

ical Imaging, 30:484–496.

BIOIMAGING 2019 - 6th International Conference on Bioimaging