FloatingCanvas : quantification of 3 D retinal structures from spectral-domain optical coherence tomography

Spectral-domain optical coherence tomography (SD-OCT) provides volumetric images of retinal structures with unprecedented detail. Accurate segmentation algorithms and feature quantification in these images, however, are needed to realize the full potential of SD-OCT. The fully automated segmentation algorithm, FloatingCanvas, serves this purpose and performs a volumetric segmentation of retinal tissue layers in three-dimensional image volume acquired around the optic nerve head without requiring any pre-processing. The reconstructed layers are analysed to extract features such as blood vessels and retinal nerve fibre layer thickness. Findings from images obtained with the RTVue-100 SD-OCT (Optovue, Fremont, CA, USA) indicate that FloatingCanvas is computationally efficient and is robust to the noise and low contrast in the images. The FloatingCanvas segmentation demonstrated good agreement with the human manual grading. The retinal nerve fibre layer thickness maps obtained with this method are clinically realistic and highly reproducible compared with time-domain StratusOCT. © 2010 Optical Society of America OCIS codes: (100.0100) Image processing; (170.4470) Ophthalmology; (170.4500) Optical coherence tomography; (170.4580) Optical diagnostics for medicine References and links 1. D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito, and J. G. Fujimoto, "Optical coherence tomography," Science 254, 1178-1181 (1991). 2. M. E. v. Velthoven, D. J. Faber, F. D. Verbraak, T. G. v. Leeuwen, and M. D. d. Smet, "Recent developments in optical coherence tomography for imaging the retina," Prog. Retin. Eye Res. 6, 57-77 (2007). 3. J. S. Schuman, M. R. Hee, A. V. Arya, T. Pedut-Kloizman, C. APuliafito, J. G. Fujimoto, and E. A. Swanson, "Optical coherence tomography: a new tool for glaucoma diagnosis," Curr Opin Ophthalmol 6, 8995 (2005). 4. M. R. Hee, J. A. Izatt, E. A. Swanson, D. Huang, J. S. Schuman, C. P. Lin, C. A. Puliafito, and J. G. Fujimoto, "Optical coherence tomography of the human retina.," Arch Ophthalmol 113, 325-332 (1995). 5. P. Carpineto, M. Ciancaglini, E. Zuppardi, G. Falconio, E. Doronzo, and L. Mastropasqua, "Reliability of nerve fiber layer thickness measurements using optical coherence tomography in normal and glaucomatous eyes," Ophthalmology 110, 190-195 (2003). 6. R. R. Bourne, F. A. Medeiros, C. Bowd, K. Jahanbakhsh, L. M. Zangwill, and R. N. Weinreb, "Comparability of Retinal Nerve Fiber Layer Thickness Measurements of Optical Coherence Tomography Instruments," Invest Ophthalmol Vis Sci 46, 1280-1285 (2005). 7. V. Guedes, J. S. Schuman, E. Hertzmark, G. Wollstein, A. Correnti, R. Mancini, D. Lederer, S. Voskanian, L. Velazquez, H. M. Pakter, T. Pedut-Kloizman, J. G. Fujimoto, and C. Mattox, "Optical Coherence Tomography Measurement of Macular and Nerve Fiber Layer Thickness in Normal and Glaucomatous Human Eyes," Ophthalmology 110, 117-189 (2003). 8. A. F. Fercher, C. K. Hitzenberger, G. Kamp, and S. Y. El-Zaiat, "Measurement of intraocular distances by backscattering spectral interferometry," Opt Commun 117, 43-48 (1995). 9. M. Wojtkowski, R. Leitgeb, A. Kowalczyk, T. Bajraszewski, and A. F. Fercher, "In vivo human retinal imaging by Fourier domain optical coherence tomography," J Biomed Opt 7, 457-463 (2002). 10. N. Nassif, B. Cense, B. H. Park, S. H. Yun, T. C. Chen, B. E. Bouma, G. J. Tearney, and J. F. d. Boer, "In vivo human retinal imaging by ultrahigh-speed spectral domain optical coherence tomography," Opt Lett 29, 480-482 (2004). 11. R. Leitgeb, C. K. Hitzenberger, and A. F. Fercher, "Performance of fourier domain vs. time domain optical coherence tomography," Opt Express 11, 889-894 (2003). 12. J. F. d. Boer, B. Cense, B. H. Park, M. C. Pierce, G. J. Tearney, and B. E. Bouma, "Improved signal-to-noise ratio in spectral-domain compared with time-domain optical coherence tomography," Opt Lett 28, 20672069 (2003). 13. H. Ishikawa, D. M. Stein, G. Wollstein, S. Beaton, J. G. Fujimoto, and J. S. Schuman, "Macular Segmentation with Optical Coherence Tomography," Invest Ophthalmol Vis Sci 46, 2012-2017 (2005). 14. D. C. Fernández, H. M. Salinas, and C. A. Puliafito, "Automated detection of retinal layer structures on optical coherence tomography images," Opt Express 13, 10200-10216 (2005). 15. C. Ahlers, C. Simader, W. Geitzenauer, G. Stock, P. Stetson, S. Dastmalchi, and U. Schmidt-Erfurth, "Automatic segmentation in three-dimensional analysis of fibrovascular pigmentepithelial detachment using high-definition optical coherence tomography," Br J Ophthalmol 92, 197-203 (2008). 16. T. Fabritius, S. Makita, M. Miura, R. Myllyl ̈a, and Y. Yasuno, "Automated segmentation of the macula by optical coherence tomography," Opt Express 17, 15659-15669 (2009). 17. M. Mujat, R. Chan, B. Cense, B. Park, C. Joo, T. Akkin, T. Chen, and J. d. Boer, "Retinal nerve fiber layer thickness map determined from optical coherence tomography images," Opt Express 13, 9480-9491 (2005). 18. D. Koozekanani, K. Boyer, and C. Roberts, "Retinal thickness measurements from optical coherence tomography using a Markov boundary model," IEEE Trans Med Imaging 20, 900-916 (2001). 19. M. Haeker, M. Sonka, R. Kardon, V. A. Shah, X. Wu, and M Abramoff, "Automated segmentation of intraretinal layers from macular optical coherence tomography images," Proc. the SPIE: Medical Imaging, 6512 (2007). 20. M. Niemeijer, M. Garvin, B. v. Ginneken, M. Sonka, and M. Abramoff, "Vessel segmentation in 3D spectral OCT scans of the retina," Proc. SPIE, 6914 (2008). 21. A. Mishra, A. Wong, K. Bizheva, and D. A. Clausi, "Intra-retinal layer segmentation in optical coherence tomography images," Opt Express 17, 23719-23728 (2009). 22. M. Garvin, M. Abramoff, R. Kardon, S. Russell, X. Wu, and M. Sonka, "Intraretinal Layer Segmentation of Macular Optical Coherence Tomography Images Using Optimal 3-D Graph Search," IEEE Trans Med Imaging 27, 1495-1505 (2008). 23. M. K. Garvin, M. D. Abramoff, X. Wu, S. R. Russell, T. L. Burns, and M. Sonka, "Automated 3-D intraretinal layer segmentation of macular spectral-domain optical coherence tomography images," IEEE Trans Med Imaging 28, 1436-1447 (2009). 24. G. Quellec, K. Lee, M. Dolejsi, M. K. Garvin, M. D. Abramoff, and M. Sonka, "Three-dimensional analysis of retinal layer texture: identification of fluid-filled regions in SD-OCT of the macula," IEEE Trans Med Imaging 29, 1321-1330 (2010). 25. K. Lee, M. Niemeijer, M. K. Garvin, Y. H. Kwon, M. Sonka, and M. D. Abramoff, "Segmentation of the optic disc in 3-D OCT scans of the optic nerve head," IEEE Trans Med Imaging 29, 159-168 (2010). 26. S. J. Chiu, X. T. Li, P. Nicholas, C. A. Toth, J. A. Izatt, and S. Farsiu, "Automatic segmentation of seven retinal layers in SDOCT images congruent with expert manual segmentation," Opt Express 18, 19413-19428 (2010). 27. E. Götzinger, M. Pircher, W. Geitzenauer, C. Ahlers, B. Baumann, S. Michels, U. Schmidt-Erfurth, and C. K. Hitzenberger, "Retinal pigment epithelium segmentation by polarization sensitive optical coherence tomography," Opt Express 16, 16410-16422 (2008). 28. K. Li, X. Wu, D. Z. Chen, and M. Sonka, "Optimal Surface Segmentation in Volumetric Images—A GraphTheoretic Approach," IEEE Trans Pattern Anal Machine Intell 28, 119-134 (2006). 29. C. K. I. Williams, "Regression with Gaussian Processes," in Mathematics of Neural Networks: Models, Algorithms and Applications, (Kluwer, 1995), 30. G. Golub and W. Kahan, "Calculating the singular value and pseudo-inverse of a matrix," SIAM Numerical Analysis 63, 205-224 (1965). 31. D. C. Hood, B. Fortune, S. N. Arthur, D. Xing, J. A. Salant, R. Ritch, and J. M. Liebmann, "Blood vessel contributions to retinal nerve fiber layer thickness profiles measured with optical coherence tomography," J Glaucoma 17, 519-528 (2008). 32. D. C. Hood, A. S. Raza, K. Y. Kay, S. F. Sandler, D. Xin, R. Ritch, and J. M. Liebmann, "A comparison of retinal nerve fiber layer (RNFL) thickness obtained with frequency and time domain optical coherence tomography (OCT)," Opt. Express 17(2009). 33. S. Jiao, R. Knighton, X. Huang, G. Gregori, and C. Puliafito, "Simultaneous acquisition of sectional and fundus ophthalmic images with spectral-domain optical coherence tomography," Opt Express 13, 444-452 (2005). 34. T. Fabritius, S. Makita, Y. Hong, R. Myllylä, and Y. Yasuno, "Automated retina shadow compensation of optical coherence tomography images," J Biomed Opt 14, 010503 (2009). 35. A. F. Frangi, W. J. Niessen, K. L. Vincken, and M. A. Viergever, "Multiscale vessel enhancement filtering," Medical Image Computing and Computer-Assisted Intervention 130-137 (1998). 36. N. V. Swindale, G. Stjepanovic, A. Chin, and F. S. Mikelberg, "Automated Analysis of Normal and Glaucomatous Optic Nerve Head Topography Images," Invest Ophthalmol Vis Sci 41, 1730-1742 (2000). 37. D. Li, D. Winfield, and D. J. Parkhurst, "Starburst: A hybrid algorithm for video-based eye tracking combining feature-based and model-based approaches," in IEEE Vision for Human-Computer Interaction Workshop at CVPR, 2005), 1-8. 38. D. L. Budenz, R. T. Chang, X. Huang, R. W. Knighton, and J. M. Tielsch, "Reproducibility of Retinal Nerve Fiber Thickness Measurements Using the Stratus OCT in Normal and Glaucomatous Eyes," Invest Ophthalmol Vis Sci 46, 2440-2443 (2004). 39. J. S. Kim, H. Ishikawa, K. R. Sung, J. Xu, G. Wollstein, R. A. Bilonick, M. L. Gabriele, L. Kagemann, J. S. Duker, J. G. Fujimoto, and J. S. Schuman, "Retinal nerve fibre layer thickness measurement reproducibility improved with spectral domain optical coherence tomography," British Journal of Ophthalmology 93:, 10571063 (2009). 40. J. C. Mwanza, R. T. Chang, D. L. Budenz, M. K. Durbin, M. G. Gendy, W. Shi, and W. J. Feuer, "Reproducibility of Peripapillary Retinal Nerve Fiber Layer Thickness and Optic Nerve Head Parameters Measured with CirrusTM HD-OCT in Glaucomatous Eyes," Invest. Ophthalmol. Vis. Sci., iovs.10-5222 (2010).


Introduction
Optical Coherence Tomography (OCT) has been widely used as a tool for evaluating the structure of the retina in cross-section [1,2]. Time-domain OCT has been used in glaucoma diagnosis and follow-up by determining retinal nerve fibre layer thickness (RNFLT) since it was adopted in the clinical setting [3][4][5][6][7]. Because of its limited speed, time-domain OCT only provides RNFLT measurements in a line scan, generally a peripapillary circle but does not provide a three-dimensional (3D) RNFLT map.
The newly developed and commercialized spectral-domain OCT (SD-OCT) [8,9] provides much faster scans [10] with improved signal-to-noise ratio [11,12] compared with time-domain OCT, for example, StratusOCT TM (Carl Zeiss Meditec, CA, USA). With this benefit, this technique represents a powerful and 'real-time' tool that potentially can be used in the clinic to assist the diagnosis and management of glaucoma. The extremely high imageacquisition speed allows a 3D image to be yielded. Each image is obtained by an in-depth axial scan (A-scan), a cross-sectional 2D scan (B-scan) consisting of a series of consecutive A-scans and the 3D image volume formed by consecutive B-scans. However, despite this high performance, such an imaging technique can only be useful clinically if there is a quantitative method to provide numerical information to the clinician. Moreover, the greater acquisition speed of SD-OCT means that a much greater amount of data is generated, which undoubtedly poses a technical challenge to the computer-assisted analysis.
Algorithms for the segmentation of images have been extensively studied since the early days of computer vision and image processing. In the studies of OCT images, some segmentation algorithms were implemented by thresholding the OCT A-scan profile, which consists of a series of 'peaks' and 'valleys' that represent various high and low tissue reflectivity [13][14][15]. A tissue layer boundary was found when the pixel intensity reached a target threshold that was sometimes adaptively chosen. These methods are computationally efficient, but are susceptible to intensity inconstancy within the individual layers. This problem was partly resolved by an efficient maximum intensity search based approach proposed by Fabritius et al [16]. In other studies, a gradient was used to form an automatic active contour to minimize the overall energy, which was defined by gradient, boundary smoothness and edge density [17]. This method is less affected by the intensity variations, but is sensitive to morphological features such as blood vessels. Koozekanani et al [18] utilized a Markov model to select and organize the edges to form a coherent boundary structure. A minimum-cost closed set approach was developed by Haeker et al [19] and Niemeijer et al [20] to identify retinal layers based on a linear combination of domain-specific cost functions. Mishra et al [21] used the image gradient to derive an external force through an adaptive kernel function and used dynamic programming to identify the continuous retinal layers within the OCT images. Garvin et al [22] attempted to find a closed set in a geometric 3-D graph that minimizes the associated costs and constraints by using an optimal graph search method (graph-cut). This method was then extended [23] such that the constraints and cost functions were learned from a training set, and by using a multiscale 3-D graph search [24,25]: the retinal surfaces were detected in a subvolume constrained by the retinal surface segmented in a down-sampled image volume. Chiu et al [26] also extended the graph-cut segmentation using dynamic programming and this helps to improve the computational efficiency when it was applied on individual B-scans. Moreover, segmentation using intrinsic tissue property such as depolarization, or polarization scrambling, on backscattered light has also been investigated using other types of OCT instruments [27].
The majority of these studies [13-17, 23, 24, 26] filtered the images to remove distractive features such as the speckle noise. These filtering methods were controlled by subjectively selected parameters and had difficulties in 'balancing' the deduction of high speckle noise and preservation of structural edges, especially in images with low contrast. Haeker et al [19] and Garvin et al [22], on the other hand, proposed image averaging to create composite images from repeat scans. The composite images had higher signal-to-noise ratio, but multiple scans (six repeat scans in this case) were needed, which may exacerbate the detrimental effects of eye movement between the scans.
A few studies performed real volumetric segmentation in 3D space [19,20,[22][23][24][25], but most other segmentation was done on individual B-scans and with the 3D layer topography created by filtering across the B-scans. The segmentation of individual B-scan is computational efficient, but the filtering may reduce sensitivity to detect a structural abnormality and the change of an abnormality over time. A key benefit of combining the information from neighbouring B-scans is that it can reduce measurement variability, especially in volumetric scans with noise or low contrast in some B-scans. Computationally speaking, the core technique behind these real 3D segmentation methods is the 3D graph search algorithm [28] that has an efficient polynomial time complexity.
The aim of this study was to develop a new segmentation algorithm, FloatingCanvas, that has a balance between the robustness and efficiency. FloatingCanvas was implemented to quantify the retinal structures in 3D image volumes around the optic nerve head (ONH) obtained with SD-OCT. It was used to process the whole image volume simultaneously and to reconstruct analytical surfaces for tissue layers or their boundaries. This method was designed to be robust to noise and artefact in image volumes and thus required no pre-processing such as filtering or image averaging. It made use of the first and higher order gradient as the natural boundary between tissue layers. In this case, the algorithm searches for the retinal pigment epithelium (RPE) and retinal nerve fibre layer (RNFL) boundaries which consequently form the RNFLT measurement. Although RNFL and RPE are in theory the two tissue layers with the strongest reflectivity in these OCT images, they and their boundaries become less identifiable in images with overall or local artefacts. FloatingCanvas was tested on images taken from both healthy and glaucomatous subjects, and was compared with manual segmentation by the human expert. It was demonstrated that the algorithm was robust enough to detect the tissue layer boundaries in images with low contrast. The RNFLT maps obtained with this method were also compared with those derived from time-domain StratusOCT TM in healthy and glaucomatous subjects.

Methods
In this study, the SD-OCT images were acquired with the RTVue-100 (Optovue, Fremont, CA, USA) using the 4mm × 4mm 3D volume scan protocol around the ONH with a depth of 2mm. This provides volumetric images with 101 B-scans comprised of 513 A-scans, each of 768 pixels in depth. Therefore, the distance between two B-scans is about 5 times that of two neighbouring A-scans.
The axis in the image is subsequently denoted as x for the direction of the B-scan, y in the direction across all B-scans and z for the direction of A-scan, and the location of a pixel in the image Im is described by a vector T , , x y z or a two-tuples ( ) , z x , where x is a column vector T , x y . The positive direction in an A-scan is defined to be from the top to the bottom of the image and will be described as 'downward' subsequently. Therefore, the value of the pixel in image Im is represented as ( ) , , Im x y z or ( ) , Im z x . The pixel coordinates were all converted to a scale in microns.

Analytical surface modelled by Gaussian Process
FloatingCanvas searches for a tissue layer, or its boundaries, in the image by deforming a 3D analytical surface that is efficiently modelled by a Gaussian Process (GP) [29]. The analytical surface is spanned by a sample of 'skeleton' points ( ) column vector containing coordinates on the x-and y-axis for the ith 'skeleton' point, and i w is the coordinate on the z-axis. The skeleton points were evenly placed along the x-and y-axis to form a regular grid in the x-y space. The interval of the 'skeleton' point grid to model the anterior and posterior RNFL boundaries was chosen to be 100µm on both x-and y-axis, and the interval for RPE was 300µm given that RPE is expected to be smoother than RNFL boundaries around ONH. The GP model acts on these 'skeleton' points to determine a function ( ) f x that provides a calculation of the surface coordinate z on z-axis for any vector coordinate x on x-and y-axis.
Similar to a Gaussian distribution, the GP is defined by a mean function and a covariance function: x is the expectation of ( ) f x , and the covariance function is defined as the expectation: In this case, the GP defines the joint probability between the skeleton points and the values * ( ) f X at arbitrary locations * X to be a Gaussian distribution: In Eq. (1), w is a column vector containing all { } 1 the columns, and a similar notation is used for * X .
where l is the length-scale of the Gaussian and is the main parameter to control the smoothness of the analytical surface, 1 x and 2 x represent the coordinates of point x on the x-and y-axis, respectively. A similar notation applies to ( , ) K X X , * ( , ) K X X and * * ( , ) K X X . 2 n δ is the prior Gaussian noise variance of w and is fixed at 100µm in the algorithm. The parameter l was set to be 150µm for the surfaces modelling the anterior and posterior RNFL boundaries, and was 450µm for the surface to detect RPE.
Using the joint probability in Eq. (1), the conditional distribution can be derived as [29]: If * X contains the coordinates of all points on the x-and y-axis, * ( ) f X returns the corresponding coordinates of the analytical surface on the z-axis.
In contrast to the conventional usage of GP where w are the fixed inputs for a regression problem, w are treated as parameters of the model in this algorithm. Therefore, the task of the algorithm is to search for the parameters w to form an analytical surface that is sufficiently close to a tissue layer or its boundaries (target surface).

Analytical surface deformation
In FloatingCanvas, an analytical surface is initialized to be at a regular location such as the top or bottom of the A-scan depending on the target surface. The analytical surface is deformed by updating the parameter w according to the forces applied on the points in the surface. Without specifying a certain target surface, the deformation process is described as a differential equation with regard to an artificial time t defined on the parameters w : where F g is a function related to a constant 'gravity' force T 0, 0, g =< > g that drives the points on the surface towards one direction along the z-axis; F Img is a force driving the surface to be attached on the target surface; Θ is a binary function consisting of one or several necessary conditions of *i x in * X being a point on the target surface, and Θ is the negation of Θ . The notation ( ) column vector, and the same notation also applies to F Img and Θ . Moreover, for notational simplicity, the resulting column vectors calculated by these three functions are simplified by removing the common variables ' ( ) Eq. (4) describes a process, as indicated by the name of the algorithm, where the analytical surface acts as a 'canvas' floating in the 3D image and is driven by different forces: the analytical surface is initially 'moved' towards the target surface by a gravity force F g ; when the analytical surface is close to the target surface (i.e. the corresponding Θ function outputs 1), the force F Img takes over F g and attaches the analytical surface on the target surface. The switch between F g and F Img is controlled by the function Θ .
To form an equation about the parameter w , the left part of Eq. (4) is expanded by inserting Eq. (2) into Eq. (4) and applying the chain rule of derivative: To form a concise notation, the matrix where † Λ is the pseudo-inverse [30] of Λ and is given by ( ) Λ acts as a projection matrix, which propagates the information from the pixels on the analytical surfaces to the skeleton points. Eq. (6) establishes the analytical surface deformation in FloatingCanvas: given the old parameter old w , the new new w can be updated by: according to the definition of the derivative in Eq. (6). The functions F g , F Img and Θ are labelled with 'old' because the input * ( ) f X in these three functions is calculated using old w , and t Δ is a sufficiently small time increment which is set to 0.1 in the algorithm.
The deformation in Eq. (7) is repeated, and in each iteration, the old w is substituted with new w in last iteration. The algorithm stops when the value of w converges, or pragmatically when the change of w becomes sufficiently small ( 0.26µm < ( 0.1 ≈ pixel on z-axis) in the implementation).

Searching for tissue layers or their boundaries
In FloatingCanvas, different tissue layers or their boundaries can be found by configuring functions F g , F Img and Θ and the 'gravity' force g . To search for the anterior boundary of RNFL, the parameter w is initialised to be 0 such that the analytical surface is at the top of the A-scan. The g is set to be T 0, 0, 30 =< > g and the functions F g , F Img and Θ are: x (10) where ' g' is the dot product; ' ∇ ' is the gradient operator when it is applied to an image. For instance, is the gradient of the image Im at the location ( ) x . This value can easily be calculated using the second derivative of the image. ∇ is the surface normal operator when it is applied to a function. For example, x , which can be analytically calculated from Eq. (2). The hat '~' above the gradient operator in ∞ (8) and (9) means that the normal has been normalised to have a length of 1.
Eq. (8) acts as a 'gravity' pulling the analytical surface downwards from the top of the Ascan. The Θ function in Eq. (10) guaranties that this gravity function stops working when the analytical surface comes across a significant gradient that is larger than g . Eq. (9) 'attaches' the surface to the local maximum of the image gradient magnitude when the function Θ outputs 1.
Similarly, the RPE layer is searched for by initialising the surface to be at the bottom of the A-scan, T 0, 0, -10 =< > g and the functions F g , F Img and Θ are: x is the maximum intensity of the pixels within a 150µm window below the surface point ( ) The force in Eq. (11) pulls the surface upwards from the bottom of the A-scan when Θ in Eq. (13) shows that the intensity on the surface is larger than the current local maximum intensity The g is set to be T 0, 0, -15 =< > g , and the functions F g and F Img are the same with Eq. (8) and (9) if *i x is not in vessel region. Otherwise, F g and F Img are set to 0. The vessel detection will be described in the subsequent section. The function Θ in Eq. (14) including the two constrains only brings the analytical surface near to the RNFL posterior boundary, which is eventually decided by the gradients in function F Img .
Benefiting from the segmentation in 3D space, the 4mm × 4mm RNFLT map can simply be calculated as the difference between the segmented anterior and posterior RNFL, without smoothing out or interpolating the individually segmented B-scans [13,14,17].
The deformation procedure illustrated in Fig. 1 demonstrates the intermediate and the final forms of the analytical surface searching for the anterior RNFL boundary. The corresponding curve in the surface was superimposed on one of the B-scans. The B-scan has been cropped so that only the region containing the tissue layers was displayed. The analytical surface is initialized as a flat panel at the top of the image volume. Before the analytical surface contacts the target boundary, it is mostly driven by the gravity g and the speckle noise in the image. Therefore, the surface becomes noisy and irregular [ Fig. 1(a), Fig. 1(d)]. It is also possible that the speckle noise and some tissues above the anterior RNFL form strong gradients that set the function Θ to be 1 before the analytical surface is close to the target boundary. An example of such an area is indicated by an arrow in Fig. 1(a). However, the deformation of the analytical surface is driven not only by the force at one point but also by the forces in the neighbouring regions. Therefore, these local exceptions that are not 'supported' by the neighbouring regions affect the deformation less and are overcome in the subsequent iterations [ Fig. 1(b)]. In Fig. 1(b), Fig. 1(e), the force ( ) the surface to the target boundary in the areas where the surface is sufficiently close to the target boundary. This procedure is repeated until all parts of the target boundary are found and the surface deformation stops [ Fig. 1(c), Fig. 1(f)]. There are four main free-parameters in FloatingCanvas: the skeleton point interval, length-scale l of the Gaussian kernel in Eq. (1), the time increment t Δ in Eq. (7) and the gravity g . The choice of these parameters represents the trade-off between the computational efficiency and the robustness of the algorithm. The skeleton point grid interval and lengthscale l are closely related and they quantify how much of the 'neighbourhood' information is taken into account for the segmentation of a particular point. The algorithm considers no information from the neighbouring area and becomes a segmentation of individual B-scans if these two parameters are very small. On the other hand, large values of these two parameters make the analytical surface too rigid to model the necessary morphology variety of the target surface. Moreover, small skeleton point interval also results in larger ( , ) K X X and † Λ matrices in Eq. (2) and (6) and the algorithm is no longer affordable by normal computing platforms. The gravity g and time increment t Δ do not affect the flexibility of the analytical surface, and they mainly control the convergence performance of the algorithm. Large values of these two parameters causes quicker convergence of the analytical surface deformation, but this risks the deformation procedure 'over-stepping' the target surface and causing convergence on the wrong surface.

Vessel detection
There has been much recent discussion about how blood vessels influence current OCT segmentation algorithms causing bias in estimates of RNFLT [31,32]. In fact, RNFLT tends to be significantly overestimated or underestimated within the area of blood vessels. It is therefore necessary to mark and delineate as far as possible the blood vessels before detecting the RNFL posterior boundary.
As it is shown in Fig. 2(b), the 'vesselness' score vn S provides a clear distinction between the vessel and non-vessel pixels, and even a simple thresholding ( 0.2 vn S > ) yields a satisfactory result. Eventually, the vessels are detected in the binary image as the area with more than 100 connected pixels.

Optic nerve head approximation
Physiologically, the RNFL near and within the ONH area changes direction and becomes the neural rim of the ONH. The idea of detecting the ONH area in this study is to exclude the ONH when displaying the RNFLT map. Therefore, the method described below was not designed to find an accurate contour line of the ONH, but to derive an approximated area of the ONH. The ONH is detected as the area bounded by the end-tips of the RPE. The RPE tip is detected using the anterior RNFL and en-face image used in vessel detection. The detected anterior RNFL surface by FloatingCanvas is approximated by the combination of a quadratic and a Gaussian surface, which is a similar to a method proposed by Swindale et al when modelling scanning laser ophthalmoscope topography [36]. The initial estimate of the ONH centre is set to be the centre of the fitted Gaussian component. The local intensity gradients at every pixel in the en-face image are then calculated in the x-y plane. The candidate pixels of RPE tips are required to meet three criteria. First, RPE tips should have a sufficiently large gradient (e.g. above the 75th percentile of the distribution of gradients). Second, considering the gradients of the intensity gradients at each pixel in the en-face image, the gradient of the intensity gradient of qualifying pixels is required to be near to 0. The gradient of the intensity gradient is near to 0 when the intensity gradient is at a local maximum. Third, there are two vectors of interest for each pixel in the en-face image: the local intensity gradient at the pixel (which is at right-angles to edges) and the vector connecting this pixel with the initial estimate of the ONH centre. It is required that the angle formed by these two vectors is smaller than 45° for the candidate RPE tip pixels. This criterion removes most pixels on vessel edges which also have strong gradient, because the angles between vectors for pixels on vessel edges are generally large (e.g. close to 90°). An ellipse is fit to the candidate RPE tip pixels using a Random Sample Consensus (RANSAC) parameter estimation (used by Li et al [37] for ellipse detection in noisy data). A cubic spine is then fitted to the pixels that are close to the fitted ellipse to form the approximated contour of the RPE tips. The ONH centre is finally calculated as the geometric mean of the contour. The ONH area is removed when the RNFLT map is displayed (Fig. 4).
The algorithms in FloatingCanvas described above were implemented under Matlab (version 7.4.0 R2007a, The MathWorks, Inc., Natick, MA).

Validation
To validate the segmentation algorithm, 26 glaucomatous subjects (mean age of 52 (range 22 to 91) years) and 14 healthy subjects (mean age of 50 (range 16 to 67) years) were recruited. The study was approved by an ethics committee and informed consent, according to the tenets of the Declaration of Helsinki, was obtained prior to examination from each subject. Each subject was imaged 3 times with the RTVue-100 system 4mm × 4mm 3D volume scan protocol. Images were acquired in both eyes of each subject during the same testing session by the same observer (PS) following the manufacturer instructions. Patient identifiers were removed from the data and the 3D volumes were transferred to a secure computer. FloatingCanvas was then applied to these 240 image volumes to extract the RNFLT measurement.
The first validation compared the automated segmentation by FloatingCanvas with segmentation by human manual grading. One of the three repeat image volumes of a randomly chosen eye from each subject was randomly selected for manual segmentation. For each selected image volume, the 101 B-scans were evenly divided into 10 sections, each of which contains about 10 B-scans. One B-scan was then randomly chosen from each section for manual segmentation. With 40 subjects, 400 B-scans were manually delineated and the segmented surface positions were compared with those produced by FloatingCanvas. The mean and standard deviation (SD) of signed and absolute difference between the manual and FloatingCanvas segmentation were then evaluated for healthy and glaucomatous subjects respectively.
The second validation hypothesis was that, if the method is reliable, the estimated RNFLT from the FloatingCanvas segmentation would be equivalent across different width annuli around the ONH. Therefore, overall mean and quadrant mean RNFLT were estimated using two different calculation annuli: wide (0.58mm wide from the inner margin radius of 1.170mm) and narrow (0.29mm wide from the inner margin radius of 1.315mm) annuli (Fig.  3). The two annuli were centred on the same circle with a radius of 1.460mm, with one annulus twice as the width of the other. Fig. 3. An illustration of the location and width of two different annuli used to calculate the mean and quadrant RNFLT. The wide (0.580mm wide from the inner margin radius of 1.170mm) annulus was twice as wide as the narrow one (0.290mm wide from the inner margin radius of 1.315mm). Both annuli were centred on the same circle with a radius of 1.460mm.
Moreover, one of the most important parameters for the quantitative analysis of imaging techniques is the reproducibility, which directly relates to the reliability of the techniques and their ability to separate physiological changes from measurement variability and also to detect progressive RNFL loss over time. The reproducibility of RNFLT measurements was evaluated by estimating test-retest variability based on the three repeated measurements and the coefficient of variation (CV) for mean and quadrant RNFLT. We defined test-retest variability of RNFLT, expressed in micrometers, as twice the SD of the three repeated measurements. The coefficient of variation was calculated as the SD of the three measurements divided by the mean.

Results
FloatingCanvas segmented the retinal structures in all 240 SD-OCT volume scans without clinically spurious results. On average, it took 5.6 1.2 ± minutes to process a large image volume ( 513 101 × A-scans) on one core of a Intel Core 2 Duo 2.4GHz CPU and 8GB RAM with single thread. Fig. 4 shows RNFLT maps from a healthy subject and a glaucomatous patient. The RNFLT map was colour-coded in micrometers with the 'hotter' colour denoting thicker RNFL. It should be noted that the healthy retina [ Fig. 4(a)] has a much thicker RNFL in the superior and inferior quadrants compared with the nasal and temporal quadrants. This is consistent with the known normal retinal anatomy. The reduced RNFLT, especially in the superior and inferior quadrants, in the glaucomatous eye [ Fig. 4(b)] can be observed and is consistent with clinical knowledge. The mean and SD of the signed and absolute difference between the manual and FloatingCanvas segmentation for glaucomatous and healthy subjects are given in Table 1. For all manually segmented B-scans, the mean and SD of the absolute difference for all three boundaries are 4.3±2.0µm. Therefore, on average, the FloatingCanvas segmentation differs from that of the human expert by 4.3µm, which is equivalent to 1.7 pixels on z-axis. The mean absolute difference between the manual and algorithm segmentation is relatively higher for glaucomatous retina compared with that of healthy retina but the difference is not statistically significant. RNFLT profile in the healthy and glaucomatous eyes is summarized in Table 2. There were no statistically significant differences between RNFLT measurements using the calculation annuli with different widths, which suggests that FloatingCanvas is robust and stable across the 3D volume. The quadrant RNFLT shows a difference between healthy and glaucomatous eyes. In general and on average, the healthy eyes, as expected, have a thicker RNFL, especially in the superior and inferior quadrants.
The reproducibility of the segmented RNFLT using SD-OCT was compared with the typical reproducibility of StratusOCT as reported in the literature using the standard scan protocol [38] in Table 3 (healthy subjects) and Table 4 (glaucomatous subjects). From Table 3 and Table 4, it can be seen that test-retest reproducibility in RNFLT measurements is better for both healthy and glaucomatous eyes with SD-OCT. RNFLT measurements were least reproducible in the nasal quadrant, with both SD-OCT and StratusOCT, while the segmented nasal measurement with SD-OCT showed markedly better reproducibility (~7µm vs 10.2µm in both normal and glaucomatous eyes). Moreover, RNFLT measurements in glaucomatous eyes were more variable than those of healthy eyes with both SD-OCT and StratusOCT, but SD-OCT showed much less variability and better reproducibility compared with StratusOCT, especially in the superior and inferior quadrants, which are the most important areas for glaucoma diagnosis. These results are consistent with the literatures about the reproducibility on another SD-OCT platform (Cirrus, Carl Zeiss Meditec, CA, USA) [39,40].

Discussion
FloatingCanvas has been developed as an effective 3D segmentation method for SD-OCT volume scans centred on the ONH. It is important that automatic segmentation should be compared with the manual segmentation as the gold standard. FloatingCanvas segmentation demonstrated good agreement with the human manual grading. It also provides a repeatable estimation of the RNFLT in the image volume. As opposed to the sparse area covered by the circular scans used in StratusOCT, the RNFLT maps cover a larger and clinically more useful area allowing for a more reliable measure of the RNFLT. The method has been tested on 240 3D volume scans acquired from both healthy and glaucomatous eyes of 40 subjects without spurious results under visual inspection. The results indicate that the RNFLT map gives a highly reproducible evaluation of a larger retina area compared with the last-generation time domain StratusOCT. The main novelty of FloatingCanvas, compared with previous algorithms, is that it processes the whole volume of data in its 'raw' format without pre-processing, such as filtering or image averaging. The algorithm needs no segmentation of individual B-scans so that it benefits from the fact that the covariance among neighbouring B-scans helps to make the analytical surface deformation robust to local noise or errors in individual B-scans. To illustrate the benefit of the volume segmentation, an individual B-scan with local low contrast is shown in Fig. 5(a). Although this image volume has good overall image quality (image quality score >60 in RTVue system), a part of the posterior RNFL boundary in this B-scan is not well defined due to the low contrast in the region marked by the white arrow in Fig. 5(a). Segmentation of the indistinct posterior RNFL boundary in this region would pose problems even for expert clinicians and B-scan-based segmentation algorithms would give spurious results in this case. However, this potential source of segmentation error can be resolved by taking into account of the neighbouring B-scans. Fig. 5(b), Fig. 5(c) show the B-scans adjacent to the one in Fig. 5(a). The locations of these three B-scans in the image volume are shown in Fig. 5(d). The adjacent B-scans have a more distinct posterior RNFL boundary in the region where it is indistinct in Fig. 5(a). The covariance among B-scans modelled by the GP model in FloatingCanvas can 'borrow' the information from the neighbouring B-scans to aid the segmentation in a local area. Therefore, the posterior RNFL boundary in Fig. 5(a) can be correctly identified by the algorithm [Fig. 5

(e)] even if the information is incomplete in this individual B-scan.
Furthermore, FloatingCanvas is computationally efficient as a segmentation algorithm in 3D space. There are two computationally intensive components in the algorithm: 1) the matrix inversion involved in the calculation of the projection matrix † Λ in Eq. (7). The time complexity of matrix inversion scales as ( ) 3 O n given the number of skeleton points n .
However, the projection matrix † Λ is not changed during the deformation of an analytical surface and thus only need to be computed once before the loop of the deformation. Practically, † Λ can be pre-computed before the algorithm and loaded into memory when it is needed because the † Λ matrix is decided only by the skeleton point interval and length-scale l which are all fixed parameters in the algorithm. This costs nearly no time on a modern computing platform; 2) the matrix multiplication and addition in Eq. (7). The right part of Eq. (7) has to be evaluated at every iteration of the deformation. However, matrix multiplication and addition are all low cost computational operations and can thus be implemented efficiently. Therefore, FloatingCanvas has a lower order polynomial computational complexity that is as efficient as the 3D graph search approach [28]. Garvin et al [23] and Quellec et al [24] have demonstrated that by using improved implementation, the execution time of the 3D graph search algorithm was significantly reduced compared with their original implementation [22]. Similarly, the current Matlab implementation of FloatingCanvas is relatively slow due to the nature of this interpreted programming language and the lack of the implementation optimisations, but the computational efficiency would be significantly improved if it is implemented using C/C++ programming language, multi-threading techniques together with the optimisation using higher order gradient information during the search.
FloatingCanvas was applied on the SD-OCT scan centred on the ONH which is designed for the assessment of the RNFLT and the optic disc; these are the standard and important examinations for the management of glaucoma. Therefore, the algorithm removed other intraretinal boundaries by applying the constraints in Eq. (14) to allow for the direct search of the posterior RNFL. The study of FloatingCanvas for segmentation of macula scans with more intra-retinal boundaries is part of our ongoing work. To search the intra-retina layers, the parameters in FloatingCanvas were altered such that the gravity force became dynamic and was decided by the gradient profile that analytical surfaces pass through during the search.
Overall, FloatingCanvas provides a robust and efficient delineation and evaluation of RNFL and RPE structures around the ONH. It can be a useful tool for clinically interpreting SD-OCT volumes for glaucoma diagnosis. The reproducible results can potentially be used for monitoring RNFLT changes in longitudinal studies. The larger scan area also improves the chance of achieving a stronger relationship of RNFLT measurements with visual function.