1.

Introduction

Optical coherence tomography (OCT) is a noninvasive imaging technique that has been used extensively on the posterior segment of the eye. The optically transparent nature of the human eye makes OCT a well-suited imaging technique for retinal imaging.¹ OCT is being increasingly used to measure the shape and thickness of the human cornea in vivo.²

Currently, ultrasound pachymetry is the primary technique used to measure the thickness of the cornea.³ Another system that can be used is the Orbscan system.⁴ Orbscan uses two slit lamps to illuminate the human eye to get information on axial curvature, elevation of the anterior and posterior surface, and corneal thickness throughout the cornea.⁵ Using OCT has two major advantages. It is a noncontact technique, meaning it is more comfortable for a patient, can be used on eyes that suffer trauma, and carries no risk of pressure on the eye altering measurements. The other major advantage is that OCT can produce images at a higher resolution and higher speed than any of the other techniques.^6,7

Anterior segment OCT (AS-OCT) allows the resolution of anterior and posterior surfaces of the entire cornea. This allows accurate measurement of the thickness and volume of the entire cornea, as well as the anterior chamber biometry, such as its angle and depth. It has several important medical applications from contact lens fitting, diagnosis and clinical evaluation, and surgical planning and monitoring to monitoring patients with eye pathologies.^8–10 In particular, obtaining accurate topography information of the anterior segment using this technique would also allow construction of patient-specific models for biomechanical modelling of the human eye.¹¹ There is currently a lack of automated measurement tools supplied with commercial OCT devices, and manual measurement is time consuming, tedious, and subject to human errors. For this reason, there is an increasing need for fully automated segmentation techniques to identify and trace anterior and posterior boundaries of the anterior segment accurately.

The segmentation of AS-OCT images has been explored in several studies. Shen et al. used a simple threshold-based model to measure the anterior surface of the cornea.¹² That study did not consider the location of the posterior boundary of the cornea. Tian et al. used a similar method to calculate the anterior chamber angle¹³ by locating the posterior boundary of the cornea near the iris. Their study did not investigate the location of the entire posterior boundary of the cornea. La Rocca et al. segmented three boundaries in the cornea using a hybrid graph theory and a dynamic programming framework.¹⁴ They were able to detect three boundaries at the center of the cornea, but their method did not segment boundaries over the entire cornea. An intelligent scissors-based method has also been used to segment five layers of the central cornea,¹⁵ but the method has two major disadvantages: It is not a fully automated method, so it still needs manual selection of initial points, and it attempts to segment only the central region of the cornea, which generally has the highest signal-to-noise ratio. To the best of our knowledge, there is no approach that can segment both the anterior and posterior boundaries of the entire cornea and the front part of the sclera in AS-OCT images.

The key challenge in segmenting the AS-OCT images is that there are regions with a low signal-to-noise ratio next to the central cornea on all images. This phenomenon, an example of which is shown in Fig. 1, is primarily due to the steepness of the cornea reducing the fringe amplitude of the signal and polarization effect, and also due to the telecentric scanning method, which reduces detection of reflected light from these regions. For this reason, the posterior surface in those areas is difficult to perceive and segment. Our observations tell us the cornea has an approximately elliptical shape. It is therefore assumed that this shape prior information can be used to address the above challenge. Previous studies have used prior shape knowledge to achieve improved segmentation in other applications^16,17 and have shown promising results on retinal OCT images. These methods used circular shape estimation to represent the shape of interest. Another previously used method for incorporating prior shape knowledge into a model is to use a good set of training images to derive the shape.¹⁸

Fig. 1

An example anterior segment OCT image acquired using the Visante AS-OCT system.

In this paper, a new technique for the automated segmentation of the anterior and posterior surfaces of the entire cornea and a small part of the sclera is presented. The technique is compared with Shen’s thresholding method¹² and Chan and Vese’s active contour without shape¹⁹ using a data set of 33 images against a reference standard built from the manual annotations made by an expert ophthalmologist (FB). These comparisons were made using three similarity measures: Dice’s coefficient, mean unsigned difference, and Hausdorff distance.²⁰

The remainder of the paper is organised as follows. Section 2 describes the dataset used in the study and the proposed segmentation technique. Section 3 presents the experimental results, and Sec. 4 discusses the results and concludes the paper.

2.

Methods

2.2.

Segmentation Framework

A three-step algorithm was developed. The first step was to preprocess the image in order to remove the central noise artefact and the iris. The next step was to obtain a coarse segmentation of the front eye using a thresholding technique. The final step used the new level set-based shape prior segmentation model to evolve the contour initialized from the coarse segmentation and achieve the final segmentation.

2.2.1.

Preprocessing step

All the AS-OCT images contain a common central noise artefact. This is intrinsic to the OCT scanning system and is caused by much higher reflection when the detector is located perpendicular to the corneal surface. This region was detected by calculating the mean intensity of each A scan of the image under consideration. The column with the highest mean intensity, and those next to it, will be considered the central noise artefact and removed by setting the intensity value of all the pixels within them to zero. The iris can complicate the shape representation essential for the final step; as such, it was detected and removed in a similar manner using the projection along the horizontal direction. The original image is shown in Fig. 2(a), while Fig. 2(b) shows the resulting image after preprocessing.

Fig. 2

Illustration of preprocessing and coarse segmentation steps. (a) Original image. (b) Preprocessed image after removing the iris and central noise artefact. (c) Coarse segmentation result.

2.2.2.

Coarse segmentation

The aim of this step is to produce an initial estimate of the corneal location (or coarse segmentation). This estimate is important because it will be used as the initial location of the curve to be evolved by the level set function in the following step, and its anterior boundary will be used to construct the shape constraint in the later stage. The technique described by Shen et al.¹² was adopted for this purpose. More specifically, an entropy filter was applied to the preprocessed image, as shown in Fig. 2(b), to produce an entropy map. The coarse segmentation is achieved by segmenting the entropy map using the Otsu’s thresholding method.²¹ Figure 2(c) shows the initial segmentation; a relatively good detection of the anterior surface can be achieved, but the posterior boundary is difficult to detect.

2.2.3.

Segmentation with level set and shape prior

In this step, the coarse segmentation will be further refined by the newly proposed level set segmentation model with shape prior. Level set techniques are widely used in image segmentation. They represent the contour of interest as the zero level set of a function valued everywhere on the image.²² This function is evolved by minimizing an energy function to achieve the segmentation.

How the energy function is constructed is important in determining the segmentation performance. In the new level set with shape prior segmentation model, the energy functional consists of a sum of three terms: the region’s fidelity term, the curve length penalty term, and the shape prior term. More specifically, the energy function is

Eq. (1)

$$E(\varphi )={\lambda }_1{E}_1(\varphi )+{\lambda }_2{E}_2(\varphi )+{\lambda }_3{E}_3(\varphi ),$$

where $E_1(\varphi )$ represents the region fidelity term, $E_2(\varphi )$ is the curve length penalty term, $E_3(\varphi )$ is the shape prior, and ${\lambda }_i$ are coefficients that determine the relative strength of each component. In particular, when ${\lambda }_3$ is 0, the equation reduces to the conventional Chan-Vese (CV) model,¹⁹ whose performance is evaluated in the next section.

For the region fidelity term $E_1(\varphi )$ in Eq. (1), an intensity-based model proposed by Chan and Vese¹⁹ has been used. The goal of the function is to split the image into two approximately homogenous regions. The region term has the form

Eq. (2)

$$E_1(\varphi )={\int }_{\mathrm{\Omega }}{[I(x,z)-u]}^{2}H(\varphi )+{[I(x,z)-v]}^2[1-H(\varphi )]\mathrm{d}x\mathrm{d}z,$$

where $I(x,z)$ is the image intensity at the pixel $(x,z)$, $u$ is the mean intensity inside the curve, $v$ is the mean intensity outside of the curve, $\mathrm{\Omega }$ is the space representing the image, and $H(\varphi )$ is the Heaviside function. The mean intensities are updated in every iteration of the model.

The curve length term, which ensures the boundary curves are smooth, favors shorter curves. A commonly used form was adopted and has the form

Eq. (3)

$$E_2(\varphi )={\int }_{\mathrm{\Omega }}\delta (\varphi )|\nabla \varphi |\mathrm{d}x\mathrm{d}z,$$

where $\delta (\varphi )$ is the regularized delta function corresponding to the gradient of Heaviside function. The shape prior term is responsible for ensuring the contour found is as close as possible to the shape prior of the cornea to be segmented. In this new model, the shape term that was incorporated into the energy function can be expressed as

Eq. (4)

$$E_3(\varphi )={\int }_{\mathrm{\Omega }}{(\varphi -{\varphi }_0)}^2\mathrm{d}x\mathrm{d}z,$$

where $\varphi$ is the level set function of the image, and ${\varphi }_0$ is the signed distance function representing the shape prior, which will be discussed below.

The location of the anterior boundary detected during the initial estimate is used to calculate the shape of the front eye. This is done by assuming the posterior boundary has a fixed relationship to the anterior boundary. First, an ellipse is fitted to the anterior boundary using a least squares fitting method,²³ and a signed distance function of the ellipse is calculated. Next, the central corneal thickness is calculated by classifying peaks in image intensity; better image quality at the center of the image means the first large peak can be assumed to be the anterior boundary, and the last peak can be assumed to be the posterior boundary. This method has been used elsewhere for central corneal thickness measurements.²⁴ Once the thickness is known, an estimate of the position of the posterior boundary can be made. The signed distance function is altered using a quadratic expression that shifts the zero point down:

Eq. (5)

$${\varphi }_{\text{lower}}(x,z)={\varphi }_{\text{upper}}(x,z)-{\varphi }_{\text{upper}}(x_t,z_t)-c_1{(x-x_t)}^2-c_2(x-x_t),$$

where ${\varphi }_0$ is the altered function, ${\varphi }_{\text{upper}}$ is the initial function based on the top surface only, $(x_t,z_t)$ is the point on the lower boundary calculated as discussed above, and $c_i$ are constants governing the strength of quadratic terms. The quadratic terms are added to account for the fact that the two surfaces of the cornea are not parallel; the posterior boundary has a greater curvature than the anterior boundary.

The product of the distance functions for the lower and upper boundaries is used to give ${\varphi }_0(x,z)$ the shape prior; that is,

Eq. (6)

$${\varphi }_0(x,z)={\varphi }_{\text{lower}}(x,z){\varphi }_{\text{upper}}(x,z),$$

where ${\varphi }_0(x,z)$ is the shape prior, while ${\varphi }_{\text{lower}}(x,z)$ and ${\varphi }_{\text{upper}}(x,z)$ are the signed distance functions corresponding to the anterior and posterior boundaries, respectively. The motivation of using ${\varphi }_0(x,z)$ is that it is a level set function representing a shape similar to the cornea. Taking the product of two signed distance functions ensures that ${\varphi }_0(x,z)$ has negative values between the boundaries and positive values everywhere else. On the boundaries, it has a value of 0. This formulation used in Eq. (3) attempts to force the level set to be sought as close as possible to the shape prior during the iterations.

2.2.4.

Minimizing the energy function

The contour is evolved toward the optimal location by minimizing the energy function described above. The Euler-Lagrange equation corresponding to the Eq. (1) was calculated. A gradient descent method was then used to solve this iteratively using the equation

Eq. (7)

$$\frac{\partial \varphi }{\partial t}=-{\lambda }_1\delta (\varphi )[{(I-u)}^2-{(I-v)}^2]+{\lambda }_2\nabla .\left(\frac{\nabla \varphi }{|\nabla \varphi |}\right)\delta (\varphi )-2{\lambda }_3(\varphi -{\varphi }_0),$$

where $t$ is an artificial time to represent the change to the level set function for each iteration. In order to speed up the analysis program, the shape constraint was updated every 20 iterations. The initial estimate described in the previous subsection was used to initialize $\varphi$.

The weighting of the different terms was determined empirically. The values used were ${\lambda }_1=1$, ${\lambda }_2=0.2$, and ${\lambda }_3=0.8$. The weighting of the terms is important, since it determines how much each particular term contributes to the overall energy function. Previous studies using level set functions have reported that changing the strength of the terms relative to the iterations produced better results.^16,17 However, the best results using this model were achieved when keeping the values fixed.

2.2.5.

Shape term with gradient (CVWSe)

In the above model, shape information is used to improve segmentation in areas where the image information alone is not enough. In an area with good image intensity, reducing the dependence on the shape may lead to improved results. To achieve this, we added to the shape term the gradient term

Eq. (8)

$${E^{\prime }}_3(\varphi )={\int }_{\mathrm{\Omega }}g{(\varphi -{\varphi }_0)}^2\mathrm{d}x\mathrm{d}z,$$

where $g$ is related to the image gradient and defined as

Eq. (9)

$$g=\frac{1}{1+\kappa |\nabla (G\otimes I)|},$$

where $\kappa$ is a constant, $(G\otimes I)$ is the convolution of the image with a Gaussian kernel to smooth the edges, and $\nabla$ is the standard Del operator that calculates the gradient of the image. When the gradient of the image is large, this decreases $g$ and results in the shape function having less effect on segmentation.

3.

Results

The four algorithms were applied to all 33 images. These were all carried out using a PC with Intel Core i5-2320 CPU @3.00 GHz and 4.00 GB RAM. The $\mathrm{mean}\pm \text{standard deviation}$ (std) time for the new algorithm was $102\pm 8\mathrm{s}$. Figure 3 shows the results from different methods overlaid on an example image, and Fig. 4 shows the segmented images achieved from different methods side by side. The $\text{mean}\pm \mathrm{std}$ of the DSC values is presented in Table 1. In particular, the mean DSC values for both of our methods are more than 0.9; this demonstrates an excellent agreement with the manual annotation. The mean value for both the CVWS and CVWSE is higher than 0.9, implying excellent agreement with the reference standard. An analysis of variance (ANOVA) test showed there is no statistically significant difference between our new models (CVWS and CVWSE; t-test, $p=0.849$). However, CVWS and CVWSE provide significantly higher DSC measures than the other two methods (CV and thresholding; all $p<0.001$).

Fig. 3

Illustration of agreement between the segmentations using the two new methods (CVWS and CVWSe) and the manual annotation. The red line is CVWS, the green line is CVWSe, and the blue line is manual annotation. Colors are altered where lines overlap. Good agreement among the different methods can be seen, especially on anterior surface.

Fig. 4

Illustration of segmentation results using different techniques: (a) CVWS, (b) CVWSe, (c) CV, (d) threshold, (e) expert, and (f) original image.

Table 1

Comparison with manual segmentation using Dice’s similarity coefficient (DSC), mean unsigned surface positioning errors (MSPE), and 95% Hausdorff distance. DSC is a coefficient; the other values are in pixels.

The MSPEs between the manual and automatic segmentation boundaries are shown in Table 1. There is no statistically significant difference in anterior boundary (ANOVA, $p=0.058$). For the posterior boundary, there is a statistically significant difference among the four methods (ANOVA, $p<0.001$). Also, there is no significant difference between CVWS and CVWSE ($p=0.212$); both methods perform significantly better than the other two methods ($p<0.001$). There is no significant difference between the CV and threshold approaches ($p=0.136$). In general, the difference is smaller for the anterior boundary, which also confirms the observation that it is much easier to detect than the posterior boundary, due to the relatively poorer image quality at the posterior cornea. For the 95% Hausdorff distance, there are statistically significant differences among the four methods ($p<0.001$). CVWSe has a larger Hausdorff distance than CVWS, but the difference is not significant ($p=0.578$); both CVWSe and CVWS have a significantly smaller Hausdorff distance than the other two methods ($p<0.001$). There are no significant differences between the CV and threshold approaches ($p=0.201$).

Our technique can be extended to segment the entire anterior chamber; the iris and remaining sclera can be easily segmented by the standard Chan and Vese model from the image where the cornea has been removed. By combining the two segmentation results together, the entire anterior chamber can be segmented. This is essential for anterior chamber biometry. Figure 5 shows an example of full segmentation, including the iris. From this figure, it appears that the anterior surface of the iris is easier to detect than the posterior surface, due to the nature of the OCT image.

Fig. 5

Example segmentation of full anterior segment, including the iris.

The ultimate goal of this work is to produce an automated technique that can detect the cornea in normal and diseased eyes. Figure 6 shows two segmentations from OCT images of an eye with keratoconus. These preliminary results show that this new segmentation technique can be used to segment the cornea in diseased eyes. It is expected that further evaluation will be performed when data from more diseased eyes are available.

Fig. 6

Two segmented images of an eye with keratoconus.

4.

Discussion and Conclusion

A fully automatic technique has been developed that can detect both the anterior and posterior surfaces of the anterior segment in AS-OCT images. The algorithm used a shape prior to allow difficult-to-segment regions to be segmented. The technique has been demonstrated to be capable of segmenting images, including regions with a low signal-to-noise ratio.

The newly developed method performed significantly better than previously described methods, and the results showed a high level of agreement with expert manual segmentation. This is the first method that has demonstrated segmentation of both anterior and posterior surfaces over the entire length of the cornea.

One of the current limitations of the algorithm is that it has not been optimized for speed. This can be done in the future by implementation in C++ or using splines to represent the level set. The work has also focused on time domain (TD)-OCT images. Although the new spectral domain (SD) OCT systems provide much faster acquisition and better resolution than TD-OCT,²⁷ there is currently no commercially available SD-OCT system that can image the entire anterior segment, including the limbus and anterior sclera. Given the advantages of SD-OCT, it is believed that the technique developed in this paper as a generic segmentation tool will be easily transferable to SD-OCT images when SD-OCT becomes mature in imaging the entire anterior segment.

One important factor that can affect segmentation performance is the image quality, including the signal-to-noise ratio (SNR). In general, the higher SNR an image has, the easier segmentation will be. For this particular problem, images contain speckle noise inherent in the OCT system and poor SNR in some of the cornea structures. This means simple thresholding and region-based models will not work; this was demonstrated in the comparative study. This study uses shape to overcome this problem. It is important to note that technical advances such as SD-OCT will lead to better-quality images. This could make segmentation easier, but this algorithm will still work easily. For SD-OCT image analysis, computational time for 3D processing will become important.

Accurate detection of the anterior and posterior surfaces is essential in research and clinical practice. For example, this information could be used as an input when creating patient-specific models of the human eye. Incorrect information or inaccurate segmentation would produce errors in the model. Patient-specific models would allow for improved diagnosis of corneal pathologies such as keratoconus and improved monitoring of the cornea after surgery. Other measurements of the cornea, such as corneal power measurements, require very precise segmentation. This could be another important application of the method developed here.

Some preliminary work presented here has shown that segmentation of patients with keratoconus is possible in principle using this technique. Future evaluation work will include investigating how well the algorithm can cope with examination of images of eyes with a variety of conditions. We also wish to extend this technique to the formation of 3D maps. This should be relatively straightforward, since the level set formulation we used extends easily to higher dimensions. In the current formulation, an intensity-based region term is used, and we plan to investigate the usefulness of texture models.^28,29

In conclusion, this work has shown that using a shape prior term can significantly improve segmentation results for the fully automatic segmentation of the cornea in AS-OCT images, a challenging and previously unsolved problem. The algorithm developed here is the first fully automatic method to detect two boundaries across the entire anterior segment with superior performance over the previous methods and excellent agreement with expert manual annotation. It may become a valuable tool for providing accurate and reliable measurements of the anterior segment geometry for clinical and nonclinical applications.