1. 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–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 image-acquisition 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–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–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 computationally 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 neighboring 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 volume images 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 artifact in volume images 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 artifacts. 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.

2. 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 neighboring 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 ${〈x,y,z〉}^T$ or a two-tuples $\left(x,z\right)$, where x is a column vector ${〈x,y〉}^T$. 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\left(x,y,z\right)$ or $Im\left(x,z\right)$. The pixel coordinates were all converted to a scale in microns.

2.1 Analytical surface modeled by Gaussian Process

FloatingCanvas searches for a tissue layer, or its boundaries, in the image by deforming a 3D analytical surface that is efficiently modeled by a Gaussian Process (GP) [29]. The analytical surface is spanned by a sample of ‘skeleton’ points ${\left\{\left(x_i\text{, }{w}_i\right)\right\}}_{i=1}^N$, where $x_i=<x_i,y_i{>}^T$ is a column vector containing coordinates on the x- and y-axis for the ith ‘skeleton’ point, and $w_i$ 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: $f(x)~GP(m(x)\text{, }k(x\text{, }{x}_{*}))$, where the mean function $m(x)=E(f(x))$ is the expectation of $f(x)$, and the covariance function is defined as the expectation: $k(x\text{, }{x}_{*})=$ $E((f(x)-m(x))(f(x_{*})-m(x_{*})))$. 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 ${\left\{w_i\right\}}_{i=1}^N$ and X is a matrix with ${\left\{x_i\right\}}_{i=1}^N$ in the columns, and a similar notation is used for $X_{*}$. $K(X\text{, }{X}_{*})$ is a matrix with element $K_{ij}=k(x_i\text{, }{x}_{*j})$, where $x_i$ and $x_{*j}$ are the ith and jth columns in X and $X_{*}$ respectively; $k(x_i\text{, }{x}_{*j})$ is defined to be an un-normalized Gaussian kernel function between data points $x_i$ and $x_{*j}$: $k(x_i\text{, }{x}_{*j})=e^{-\frac{1}{2}\left(\frac{{(\left(x_i^1-x_{*j}^1\right))}^2+{\left(x_i^2-x_{*j}^2\right)}^2}{l^2}\right)}$, where l is the length-scale of the Gaussian and is the main parameter to control the smoothness of the analytical surface, $x^1$ and $x^2$ represent the coordinates of point x on the x- and y-axis, respectively. A similar notation applies to $K(X\text{, }X)$, $K(X_{*}\text{, }X)$ and $K(X_{*}\text{, }{X}_{*})$. ${\delta }_n^2$ 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 $p(f(X_{*})|X\text{, }{X}_{*}\text{, }w)$ $~N(\overline{f}(X_{*})\text{, cov(}\overline{f}(X_{*})\text{)}$ can be derived as [29]:

If $X_{*}$ contains the coordinates of all points on the x- and y-axis, $\overline{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 input 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 parameter w to form an analytical surface that is sufficiently close to a tissue layer or its boundaries (target surface).

2.2 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:

Equation (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 ${\text{F}}_g$; when the analytical surface is close to the target surface (i.e. the corresponding Θ function outputs 1), the force ${\text{F}}_{Img}$ takes over ${\text{F}}_g$ and attaches the analytical surface on the target surface. The switch between ${\text{F}}_g$ and ${\text{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 $K(X_{*}\text{, }X){(K(X\text{, }X)+{\delta }_n^{2}I)}^{-1}$ is substituted by Λ. Multiplying ${\left({\mathrm{\Lambda }}^{\text{T}}\mathrm{\Lambda }\right)}^{-1}{\mathrm{\Lambda }}^{\text{T}}$on both side of Eq. (5) gives:

Equation (6) establishes the analytical surface deformation in FloatingCanvas: given the old parameter $w^{old}$, the new $w^{new}$ can be updated by:

The deformation in Eq. (7) is repeated, and in each iteration, the $w^{old}$ is substituted with $w^{new}$ in last iteration. The algorithm stops when the value of w converges, or pragmatically when the change of w becomes sufficiently small ($<0.26\text{μm}$ ($\approx 0.1$ pixel on z-axis) in the implementation).

2.3 Searching for tissue layers or their boundaries

In FloatingCanvas, different tissue layers or their boundaries can be found by configuring functions ${\text{F}}_g$, ${\text{F}}_{Img}$ and Θ and the ‘gravity’ force g. To search for the anterior boundary of RNFL, the parameter w is initialized to be 0 such that the analytical surface is at the top of the A-scan. The g is set to be $g=<0,0,\text{ 30}{>}^{\text{T}}$ and the functions ${\text{F}}_g$, ${\text{F}}_{Img}$ and Θ are:

Equation (8) acts as a ‘gravity’ pulling the analytical surface downwards from the top of the A-scan. 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 $\left|g\right|$. Equation (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 initializing the surface to be at the bottom of the A-scan, $g=<0,0,\text{ -10}{>}^{\text{T}}$ and the functions ${\text{F}}_g$, ${\text{F}}_{Img}$ and Θ are:

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 $I{m}_{\max }\left(x_{*i}\right)$. Equation (12) ‘attaches’ the surface to the locations with the maximum local intensity when the function Θ outputs 1.

To search for the posterior RNFL, the surface is initialized to be the detected RPE. Two necessary condition constraints are applied during the search of the posterior RNFL: the intensity on the posterior RNFL should be no less than 70% of the intensity of the RNFL anterior boundary, and the gradient of the RNFL posterior boundary should be ‘downward’ (${\nabla }_{z}Im\left(x_{*i},\overline{f}(x_{*i})\right)$>0). These two constraints are integrated into the function Θ:

The g is set to be $g=<0,0,\text{ -15}{>}^{\text{T}}$, and the functions ${\text{F}}_g$ and ${\text{F}}_{Img}$ are the same with Eq. (8) and (9) if $x_{*i}$ is not in vessel region. Otherwise, ${\text{F}}_g$ and ${\text{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 ${\text{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), 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 neighboring regions. Therefore, these local exceptions that are not ‘supported’ by the neighboring regions affect the deformation less and are overcome in the subsequent iterations [Fig. 1(b)]. In Fig. 1(b), 1(e), the force ${\text{F}}_{img}\left(x_{*i},\overline{f}(x_{*i})\right)$ from Eq. (9) ‘attaches’ 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), 1(f)].

 

Fig. 1 Steps in the deformation procedure to search for the anterior RNFL boundary. (a-c) a B-scan in an image volume superimposed with the corresponding curve in the analytical surface at the 30th, 45th and the last (96th) iterations of the deformation procedure. Plates (d-f) show an overview of the analytical surface at these 3 iterations. The arrow in (a) shows an area with strong gradient that sets the Θ function to be 1 before the analytical surface is close to the target boundary.

Download Full Size | PDF

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 $\mathrm{\Delta }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 length-scale l are closely related and they quantify how much of the ‘neighborhood’ information is taken into account for the segmentation of a particular point. The algorithm considers no information from the neighboring 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\text{, }X)$ and ${\mathrm{\Lambda }}^{\text{†}}$ matrices in Eq. (2) and (6) and the algorithm is no longer affordable by normal computing platforms. The gravity g and time increment $\mathrm{\Delta }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.

2.4 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.

FloatingCanvas identifies blood vessels by using the en-face image $EF(x)$ obtained by averaging the 50 pixels below and above the analytical RPE [33,34] [Fig. 2(a) ]. This detection scheme computes the ‘vesselness’ for each pixel in the en-face image. ‘Vesselness’ is a definition based on the analysis of eigenvalues of the Hessian matrix of image intensity [35]:

 

Fig. 2 En-face image and pixel vesselness. (a) the en-face image calculated by averaging the 50 pixels below and above the detected RPE. (b) the pixel vesselness in grayscale.

Download Full Size | PDF

2.5 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 ).

 

Fig. 4 The RNFLT map calculated from the segmented retinal structures of a healthy subject (a) and a glaucomatous patient (b). The RNFLT map was colour-coded in micrometers. Note the significantly thicker RNFLT, especially in the superior and inferior areas, in the healthy example (a) compared with that of the glaucomatous example (b).

Download Full Size | PDF

The algorithms in FloatingCanvas described above were implemented under Matlab (version 7.4.0 R2007a, The MathWorks, Inc., Natick, MA).

2.6 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 volume images 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 was 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 centered 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 centered on the same circle with a radius of 1.460mm.

Download Full Size | PDF

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.

3. Results

FloatingCanvas segmented the retinal structures in all 240 SD-OCT volume scans without clinically spurious results. On average, it took $5.6\pm 1.2$ minutes to process a large image volume ($513\times 101$ A-scans) on one core of a Intel Core 2 Duo 2.4GHz CPU and 8GB RAM with single thread.

Figure 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 .

Table 1. Mean and SD of signed and absolute difference between the manual and FloatingCanvas segmentation for glaucomatous and healthy subjects. The difference values were summarized with all 3 surfaces and anterior, posterior RNFL and RPE respectively

View Table | View all tables in this article

Table 2. Mean and SD of total and quadrant retinal nerve fiber layer thickness (RNFLT) of healthy and glaucomatous eyes. RNFLT was determined with two types of calculation annuli (0.58mm and 0.29mm wide, respectively)

View Table | View all tables in this article

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].

Table 3. Coefficient of variation and test-retest variability of total and quadrant retinal nerve fiber layer thickness (RNFLT) of healthy retina measured by FloatingCanvas. RNFLT was calculated at two widths of calculation annuli (0.58mm and 0.29mm). A typical reproducibility of StratusOCT is given for comparison

View Table | View all tables in this article

Table 4. Coefficient of variation and test-retest variability of total and quadrant retinal nerve fiber layer thickness (RNFLT) of glaucomatous subjects measured by FloatingCanvas. RNFLT was calculated at two widths of calculation annuli (0.58mm and 0.29mm). A typical reproducibility of StratusOCT is given for comparison

View Table | View all tables in this article

4. Discussion

FloatingCanvas has been developed as an effective 3D segmentation method for SD-OCT volume scans centered on the ONH. It is important that automatic segmentation should be compared with the manual segmentation as the gold standard. The 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 neighboring 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 neighboring B-scans. Figure 5(b), 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 modeled by the GP model in FloatingCanvas can ‘borrow’ the information from the neighboring 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.

 

Fig. 5 The benefit of a segmentation in 3D space. (a) a B-scan with a segment of RNFL with an indistinct posterior boundary. The region with an indistinct posterior RNFL boundary is denoted by a white arrow in the image. (b, c) the RNFL boundary in the B-scans adjacent to the one in (a) is more distinct; the location of indistinct RNFL boundary in (a) is also marked. (d) the locations of B-scans (a, b, c) in the image volume. (e) the segmented RPE and anterior and posterior RNFL boundaries in image (a) by FloatingCanvas.

Download Full Size | PDF

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 ${\mathrm{\Lambda }}^{\text{†}}$ in Eq. (7). The time complexity of matrix inversion scales as $O\left(n^3\right)$ given the number of skeleton points n. However, the projection matrix ${\mathrm{\Lambda }}^{\text{†}}$ 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, ${\mathrm{\Lambda }}^{\text{†}}$ can be pre-computed before the algorithm and loaded into memory when it is needed because the ${\mathrm{\Lambda }}^{\text{†}}$ 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 optimizations, 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 centered 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 intra-retinal 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.