Fast segmentation of anterior segment optical coherence tomography images using graph cut
© Williams et al.; licensee BioMed Central. 2016
Received: 26 August 2014
Accepted: 11 January 2015
Published: 22 January 2015
Optical coherence tomography (OCT) is a non-invasive imaging system that can be used to obtain images of the anterior segment. Automatic segmentation of these images will enable them to be used to construct patient specific biomechanical models of the human eye. These models could be used to help with treatment planning and diagnosis of patients.
A novel graph cut technique using regional and shape terms was developed. It was evaluated by segmenting 39 OCT images of the anterior segment. The results of this were compared with manual segmentation and a previously reported level set segmentation technique. Three different comparison techniques were used: Dice’s similarity coefficient (DSC), mean unsigned surface positioning error (MSPE), and 95% Hausdorff distance (HD). A paired t-test was used to compare the results of different segmentation techniques.
When comparison with manual segmentation was performed, a mean DSC value of 0.943 ± 0.020 was achieved, outperforming other previously published techniques. A substantial reduction in processing time was also achieved using this method.
We have developed a new segmentation technique that is both fast and accurate. This has the potential to be used to aid diagnostics and treatment planning.
KeywordsOptical coherence tomography Anterior segment Image segmentation Cornea Graph cut Shape prior
The human eye is a remarkable pressurised organ and its biomechanical properties are essential in maintaining its functionality. There has been an increasing interest in modelling the biomechanics of the eye, as this will improve our understanding and management of eye disease. When carrying out biomechanical modelling of the eye, one of the aims is to be able to produce patient specific models so that the eyes of individual patients can be simulated providing important information towards personalised medicine. In order to produce patient specific models, geometry of the eyes in particular the cornea is required. A technique that can be used to gain such information on the cornea is optical coherence tomography (OCT) imaging. OCT is widely used in retinal imaging  and has increasing uses in imaging the anterior segment .
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, surgical planning and monitoring, to monitoring patients with eye pathologies [3-5]. 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 . A potential use of biomechanical modelling is for studying patients with keratoconus [7,8]. Keratoconus causes a deformation of the cornea due to change in stiffness of parts of the cornea. Currently, surgical treatment is limited to making changes to the entire cornea . A more targeted patient specific intervention in this disease has the potential to improve the treatment.
There has been some previous preliminary work in using OCT images as an input for modelling. One group studying acute angle-closure glaucoma (AACG) made measurements using OCT images and then used these to create a finite element model of the anterior chamber . Their work managed to successfully create a model, which could give information about patients suffering from AACG. This shows the potential use for AS-OCT images in biomechanical modelling. However, their underlying program for the extraction of the structures of the anterior segment is not evaluated for accuracy.
The automated segmentation of the cornea in OCT images is a prerequisite step to derive the geometry of the cornea for biomechanical modelling. Of the currently commercially available OCT devices, only CASIA (Tomey, Tokyo, Japan) is able to both image and automatically segment the entire anterior segment . However, we do not know the proprietary details of the segmentation approach used. The use of proprietary software only coming with the device will also limit the access for many applications. Other devices are either limited to only imaging a small section of the cornea or are unable to automatically segment the entire anterior segment like the Visante OCT system (Zeiss Meditec Inc., Dublin, California) . The widespread use of OCT for biomechanical modelling requires imaging to be carried out by whatever device is present. Current automated measurement tools supplied with Visante OCT devices are limited to the central region of the cornea only, and manual measurement is time consuming, tedious and subject to human errors. For this reason, the development of fully automated segmentation techniques to accurately identify and trace both anterior and posterior boundaries of the anterior segment is required.
Segmentation of AS-OCT images has previously been carried out using a threshold-based technique , but this method was unable to accurately detect the posterior surface of the cornea. Work has been recently carried out using level set segmentation techniques on OCT images of the human cornea . This technique achieved segmentation results comparable to those carried out by a manual observer. This method however, suffers from problems of relatively slow speed. Another group used graph theory and dynamic programming to segment a set of spectral domain OCT (SD-OCT) images of the central section of the cornea . This method used the image gradient to identify edges on the images; Dijkstra’s algorithm  was then used to find a shortest path through the image based on the image gradient. Their results had good agreement with manual observers only for the central region of the cornea where the highest signal to noise ratio was found. Regions further away from the centre were approximated by fitting a polynomial curve to extend the surface. This method would not be very reliable if accurate information of a wider area of the cornea is wanted.
We propose a graph cut based segmentation technique for efficient and accurate segmentation of the cornea. Graph cuts are a way of segmenting an image exploiting a min-cut algorithm for graph labelling . Graph cut methods have been used increasingly in image segmentation and other image processing problems. A set of pixels and a set of labels are input to the model and the goal is to minimize an energy function. A global optimum cut can be obtained for an image. Other groups have used graph cut segmentation to segment OCT images of the human retina [18,19]. The motivation of using this technique is that there are very efficient max flow algorithms that can be used to minimize the above energy . In order to use graph cuts, an energy function must be constructed in the correct form . There are many different terms that can be included into this function, such as terms based on regional statistics , and knowledge of shape .
In this paper, a new technique for segmentation of anterior and posterior surfaces of AS-OCT images is presented. Graph cut based segmentation is used to improve the speed and efficiency of the segmentation technique when compared to previously reported segmentation techniques. The results of the segmentation are evaluated over a data set, which has been segmented manually by an expert ophthalmologist.
39 AS-OCT B scan images through the centre of the cornea from healthy eyes (one image per eye) were used for the study. The study was approved by the Institutional Review Board and undertaken by following the tenets of the Helsinki Declaration. These images were acquired using the Visante AS-OCT system in Wenzhou Medical University, China. The Visante system is a time domain system that uses 1,300 nm infrared light to obtain cross-sectional images of the anterior segment with a scanning rate of 2,000 axial scans per second. Each B scan image contains 256 A-scans in 16 mm with 1024 points per A scan to a depth of 8 mm. The images have a transverse resolution of 60 μm and an axial resolution of 18 μm. The images were manually delineated by an expert ophthalmologist (FB) who marked the locations of the anterior and posterior boundaries.
Step 1: Removing iris and central noise
We have two structures present shown in our images; the cornea that is a curved structure and the iris that is largely horizontal in nature and has a gap in the centre. Both the iris and the cornea are defined by an increased intensity in the OCT images due to their dense structure scattering more light. This means that the iris has the potential to interfere with the detection of the cornea, so the exclusion of the iris from the image is desirable. OCT works by detecting scattered light. Since reflected light from the centre is of a much higher intensity than scattered light from within the cornea, most of OCT images will contain a bright, superior-inferior, narrow region at their centres. This region is an artefact due to the specula reflection of light from the centre where the cornea is perpendicular to the light source.
Both the iris and central noise can be detected and removed using a similar method. The central noise artefact is found by looking at the mean intensity of each column of pixels in the image. The columns with the maximum mean intensity correspond to the area of the central noise. Once this is detected, the intensities of the pixels within this region are all set to zero and excluded from the subsequent analysis. The same technique can be used for the detection of the iris. This time, the mean intensity of the rows is found. All rows below the horizontal line that is 10 pixels above the maximum intensity row are cropped out of the image. This line was empirically chosen since it resulted in the iris being cropped off in most images.
Step 2: Initial segmentation using a threshold
The purpose of this step is to generate an approximate shape that is used to help guide the graph cut segmentation step. First, an entropy filter was applied to the image; this gave a measure of the texture of the image. The entropy image increased the contrast between the cornea and the background allowing a threshold to be used that would find most of the cornea. The value of the threshold was chosen using Otsu’s technique , which attempts to maximize variance between the two different classes while minimising the variance within each class at the same time. Once the threshold was chosen, all pixels with a value higher than the threshold were labelled as the cornea and the rest as part of the background.
An ellipse was then fitted to the anterior surface of the segmented shape. This was done since the cornea shape can be approximated as an ellipsoid. A related ellipse was then generated to give an initial estimate of the posterior surface since the threshold method could detect the anterior surface much better than the posterior surface. This second ellipse was generated by shrinking the first ellipse so that its top edge shifted by the distance of the central cornea width, which was calculated using an existing technique . These two ellipses were then combined to generate an approximate shape for the cornea, which is used to help guide the next step of the segmentation.
Step 3: Graph Cut Segmentation
The novel graph cut image segmentation program is used to segment the images. To use this to segment an image, we constructed a graph from the image. Each pixel on the image corresponds to a point on the graph. The values of the different points on the graph are controlled  to make the optimum path correspond to one that separates the cornea from the background to achieve segmentation.
Our segmentation method creates a function with three different terms to control the weighting of the points on the graph. The three terms used were a regional statistics term, a curvature term, and a shape term. The regional statistics term will act to segment the image according to the brightness of the pixels. It will find a segmentation that splits the light pixels from the dark pixels. Acting alone, this term will produce results similar to the threshold technique used in the initial step. The curvature term acts to keep the boundary smooth. It favours shorter boundaries by smoothing out any small boundary deviations due to image noise. The shape term is generated using the shape created in the second step to encourage the segmentation to be close to the boundaries of the prior shape. This helps the detection of parts of the cornea with low signal to noise ratio. Once the energy function was constructed, the minimum cut of this function could be found very quickly using the Dynamic Boykov-Kolmogorov algorithm .
When testing the segmentation program, the strength of the regional and curvature terms were kept fixed. The optimum strength of the shape term was empirically found using a leave-one-out test for each image.
A comparison was carried out between the results of segmentation using this method and manual segmentation. Three different measures were used to evaluate the segmentation Dice’s similarity coefficient (DSC), mean unsigned surface positioning error (MSPE), and the Hausdorff distance (HD). DSC measures the overlap of two areas and gives a value of 1 for perfect agreement . MSPE measures the unsigned distance between the boundaries from the program and manual segmentation. 95% HD measures the maximum distance between each point in one boundary and the corresponding point in the other. 95% HD excludes the longest 5% of the distances before finding the maximum distance between points in the reference and actual data .
A comparison of the results of different strengths of the shape constraint was carried out using the DSC. A comparison was also made with the results of a previous published level set based technique . All segmentation was carried out using a Win7 PC with Intel Core i5-2320 CPU @ 3.00 GHz and 4.00 GB RAM. In order to test the significance of any differences, a paired t-test was carried out using SPSS (version 20, IBM). A p value of <0.05 is considered statistically significant.
Dice’s similarity coefficient (DSC) comparison between segmentation results of graph cuts program with the manual segmentation over all the images
Value of shape (weighting term)
Comparison between segmentation results and manual segmentation
MSPE - posterior
Mean time taken for segmentation using graph cut method (GC) and level set method (LS)
A new segmentation technique for anterior segmentation OCT images has been developed. The technique has demonstrated high-speed segmentation with an average segmentation time of 2.53 seconds per image. The segmentation technique has shown comparable accuracies to those achieved using manual segmentation in a fraction of the time it would take to manually segment an image.
Our new technique was compared to a previously published level set based technique. Using four different comparison methods, the newly developed method was shown to be significantly better. Another major advantage of our new method is the speed of the technique. The technique presented here is over 50 times faster. This difference is very important when looking at practical applications of the technology. If the technique is to be useful as a diagnostic tool, it is essential to process the image fast enough, so that it can be completed while the patient is still present.
Accurate segmentation of the cornea is important for analysis of structures further into the eye. Light will experience refraction when it passing through the cornea, which will affect the detection of other structures in the image. Since our focus of segmentation was to investigate the cornea, we have not carried out any correction of the image. If accurate information on the iris was required, then this method could be used to calculate the distortion that the cornea will cause to the rest of the image, allowing a correction to be made.
Further improvements to the model could be made by carrying out further investigations into the strength of the shape term to use over a larger data set. A different method of representing the shape could produce improved results for this technique, such as the use of a statistical shape model  created using a training set of manually segmented images. While this method has shown significant increases in speed compared to previous techniques, it might be possible to further increase the speed by a process of code optimisation. This segmentation has been carried out using a combination of Matlab and C++. Exclusive use of C++ may enable further speedup.
Another improvement that could be made is to alter the step involving removal of the iris. Currently, the step works by removing an entire section of the image where the iris is present. In cases where the iris is deformed, such as by iris tumours or other problems, this may result in too large a section of the image being removed. One solution to this would be to use a more advanced iris removal technique. For example, we have previously demonstrated that segmentation of the iris is much easier than segmentation of the cornea . Segmentation of the iris could be used to remove only the iris from the part of the image, which it is present in.
Future work on this technique will focus on expanding the segmentation technique to full 3D segmentation. This will then enable our segmentation to be used to create patient specific biomechanical models of the cornea.
This is the first time graph cut segmentation has been demonstrated for use on AS-OCT images. The method has been demonstrated to be faster than previous techniques with improved accuracy. The speed of this technique means it can be used to segment images while a patient waits and allow for rapid interpretation of the images. This technique can thus be used to produce an input for patient specific biomechanical models of the human eye.
Appendix: segmentation energy function
Where λ image , and λ image are coefficients determining the strength of image and shape terms respectively, u i is image intensity at point i, c1 is mean intensity of object, c2 mean intensity of background, x i is binary function that determines if point i is part of the object or background and M (p) is a signed distance function defined by the shape that guides the segmentation.
- Swanson EA, Izatt JA, Hee MR, Huang D, Lin CP, Schuman JS, et al. In vivo retinal imaging by optical coherence tomography. Opt Lett. 1993;18(21):1864–6.View ArticlePubMedGoogle Scholar
- Kaluzy BJ, Kałuzny JJ, Szkulmowska A, Gorczyńska I, Szkulmowski M, Bajraszewski T, et al. Spectral optical coherence tomography: A novel technique for cornea imaging. Cornea. 2006;25(8):960–5.View ArticlePubMedGoogle Scholar
- Sakata LM, Lavanya R, Friedman DS, Aung HT, Gao H, Kumar RS, et al. Comparison of gonioscopy and anterior segment ocular coherence tomography in detecting angle closure in different quadrants of the anterior chamber angle. Ophthalmology. 2008;115(5):769–74.View ArticlePubMedGoogle Scholar
- Konstantopoulos A, Kuo J, Anderson D, Hossain P. Assessment of the use of anterior segment optical coherence tomography in microbial keratitis. Am J Ophthalmol. 2008;146(4):534–42. e2.View ArticlePubMedGoogle Scholar
- Hall RC, Mohamed FK, Htoon HM, Tan DT, Mehta JS. Laser in situ keratomileusis flap measurements: Comparison between observers and between spectral-domain and time-domain anterior segment optical coherence tomography. J Cataract Refract Surg. 2011;37(3):544–51.View ArticlePubMedGoogle Scholar
- Elsheikh A, Wang D. Numerical modelling of corneal biomechanical behaviour. Comput Methods Biomech Biomed Engin. 2007;10(2):85–95.View ArticlePubMedGoogle Scholar
- Gefen A, Shalom R, Elad D, Mandel Y. Biomechanical analysis of the keratoconic cornea. J Mech Behav Biomed Mater. 2009;2(3):224–36.View ArticlePubMedGoogle Scholar
- Carvalho LA, Prado M, Cunha RH, Neto AC, Paranhos Jr A, Schor P, et al. Keratoconus prediction using a finite element model of the cornea with local biomechanical properties. Arq Bras Oftalmol. 2009;72(2):139–45.View ArticlePubMedGoogle Scholar
- Raiskup-Wolf F, Hoyer A, Spoerl E, Pillunat LE. Collagen crosslinking with riboflavin and ultraviolet-A light in keratoconus: Long-term results. J Cataract Refract Surg. 2008;34(5):796–801.View ArticlePubMedGoogle Scholar
- Su L, Ju Y, Liu X. Quantitative modeling and simulation of Anterior Chamber in OCT images. In: Proceedings 2010 IEEE 5th International Conference on Bio-Inspired Computing: Theories and Applications, BIC-TA 2010. 2010. p. 1314–18.Google Scholar
- Liu S, Yu M, Ye C, Lam DSC, Leung CKS. Anterior chamber angle imaging with swept-source optical coherence tomography: An investigation on variability of angle measurement. Invest Ophthalmol Vis Sci. 2011;52(12):8598–603.View ArticlePubMedGoogle Scholar
- Reinstein DZ, Gobbe M, Archer TJ. Anterior segment biometry: A study and review of resolution and repeatability data. J Refract Surg. 2012;28(7):509–20.View ArticlePubMedGoogle Scholar
- Shen M, Cui L, Li M, Zhu D, Wang MR, Wang J. Extended scan depth optical coherence tomography for evaluating ocular surface shape. J Biomed Opt. 2011;16(5):056007. doi:10.1117/1.3578461.PubMed CentralView ArticlePubMedGoogle Scholar
- Williams D, Zheng Y, Bao F, Elsheikh A. Automatic segmentation of anterior segment optical coherence tomography images. J Biomed Opt. 2013;18(5):056003. doi:10.1117/1.jbo.18.5.View ArticleGoogle Scholar
- LaRocca F, Chiu SJ, McNabb RP, Kuo AN, Izatt JA, Farsiu S. Robust automatic segmentation of corneal layer boundaries in SDOCT images using graph theory and dynamic programming. Biomed Opt Express. 2011;2(6):1524–38.PubMed CentralView ArticlePubMedGoogle Scholar
- Dijkstra EW. A note on two problems in connexion with graphs. Numerische Mathematik. 1959;1(1):269–71.View ArticleGoogle Scholar
- Boykov Y, Veksler O, Zabih R. Fast approximate energy minimization via graph cuts. IEEE Trans Pattern Anal Mach Intell. 2001;23(11):1222–39.View ArticleGoogle Scholar
- Haeker M, Sonka M, Kardon R, Shah VA, Wu X, Abràmoff MD. Automated segmentation of intraretinal layers from macular optical coherence tomography images. Progress in Biomedical Optics and Imaging - Proceedings of SPIE. 2007.Google Scholar
- Song Q, Bai J, Garvin MK, Sonka M, Buatti JM, Wu X. Optimal multiple surface segmentation with shape and context priors. IEEE Trans Med Imaging. 2013;32(2):376–86.PubMed CentralView ArticlePubMedGoogle Scholar
- Boykov Y, Kolmogorov V. An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision. IEEE Trans Pattern Anal Mach Intell. 2004;26(9):1124–37.View ArticlePubMedGoogle Scholar
- Kolmogorov V, Zabih R. What Energy Functions Can Be Minimized via Graph Cuts? IEEE Trans Pattern Anal Mach Intell. 2004;26(2):147–59.View ArticlePubMedGoogle Scholar
- El-Zehiry N, Xu S, Sahoo P, Elmaghraby A. Graph cut optimization for the mumford-shah model. In: 7th IASTED International Conference on Visualization, Imaging, and Image Processing; 2007. Spain: Palma de Mallorca; 2007. p. 182–7.Google Scholar
- Slabaugh G, Unal G. Graph cuts segmentation using an elliptical shape prior. In: IEEE International Conference on Image Processing. Italy: Genova; 2005. p. 1222–5.Google Scholar
- Otsu N. A Threshold selection method from gray-level histograms. IEEE Trans Syst Man Cybern. 1979;SMC-9(1):62–6.Google Scholar
- Bechmann M, Thiel MJ, Roesen B, Ullrich S, Ulbig MW, Ludwig K. Central corneal thickness determined with optical coherence tomography in various types of glaucoma. Br J Ophthalmol. 2000;84(11):1233–7.PubMed CentralView ArticlePubMedGoogle Scholar
- Daněk O, Matula P, Maška M, Kozubek M. Smooth Chan-Vese segmentation via graph cuts. Pattern Recognit Lett. 2012;33(10):1405–10.View ArticleGoogle Scholar
- Zou KH, Warfield SK, Bharatha A, Tempany CMC, Kaus MR, Haker SJ, et al. Statistical Validation of Image Segmentation Quality Based on a Spatial Overlap Index. Acad Radiol. 2004;11(2):178–89.PubMed CentralView ArticlePubMedGoogle Scholar
- Huttenlocher DP, Klanderman GA, Rucklidge WJ. Comparing images using the Hausdorff distance. IEEE Trans Pattern Anal Mach Intell. 1993;15(9):850–63.View ArticleGoogle Scholar
- Bresson X, Vandergheynst P, Thiran JP. A variational model for object segmentation using boundary information and shape prior driven by the Mumford-Shah functional. Int J Comput Vis. 2006;68(2):145–62.View ArticleGoogle Scholar
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.