This study was approved by the ethics committee of the Cixi Institute of Biomedical Engineering, Chinese Academy of Sciences, and adhered to the principles of the Declaration of Helsinki. Written informed consent was obtained from each subject before they participated in the study.

Two datasets were used in the experiments. Dataset1 consists of 37 anterior segment OCT (AS-OCT) B-scan images of healthy eyes acquired with a Visante AS-OCT system [see Williams et al. (13) for details]. Briefly, each image covers a 16 mm wide region sampled by 256 A-scans of 1024 points to a depth of 8 mm. The pixel resolution is therefore 60μm×18μm. All images were manually delineated by two ophthalmologists, one of whom marked the images twice in a masked fashion. All 37 images were used as test data for validation.

Dataset 2 is a 3D volume scan of a porcine cornea acquired with an in-house spectral domain OCT device at a scan rate of 100 μs per A-scan using a light source with a central wavelength of 840 nm. It consists of 421 raster B-scan images of a 15.1 mm region. The images have an axial resolution of 1.9 μm and a transverse resolution of 15 μm with a gap of 20 μm between consecutive B-scans. Three visible layer interfaces (air-epithelium, epithelium-stroma and endothelium-aqueous) of randomly selected 6 B-scans were manually marked twice by an experienced grader for validation purposes.

All images in dataset 1 generated by AS-OCT were used to evaluate the performance of the algorithm. These images were acquired using the Visante AS-OCT system, which is a time-domain system that acquires images at 2,000 axial scans per second in 1,300 nm infrared light. Each B-scan, 37 B-scans in total, consists of 816 A-scans and 406 lateral pixels in each A-scan. The scan width and depth are 16 mm and 8 mm respectively. The pixel resolution is therefore 19.70 μm × 19.60 μm. The layer interfaces for all 37 test datasets were marked by three different graders simultaneously. Test data for the second dataset, obtained from a home-made OCT, was generated by randomly selecting 6 B-scans from a pool of 421 OCT images. The layer boundaries for the test data were then delineated by an experienced grader. The scan rate of the home-built spectral domain OCT system is 100 ms in each A-scan using a light source with a central wavelength of 840 nm. Each B-scan image contains 1000 A-scans of 15.11 mm with 1024 pixels of 1.9511 mm. The pixel resolution for the second dataset is therefore 1.905 μm × 15.113 μm.

The first step is to remove unwanted structures and noise (e.g., the iris and high intensity artifacts at the apex position) from the image, as shown in Figure 2, to reduce their detrimental effect on segmentation performance. This is achieved by cropping the image, after which all the content remaining in the image becomes a region of interest (ROI). It is observed that the apex of the cornea in OCT images has a relatively higher intensity than the region above it, as the scattered light from other regions is much weaker than that in the center. Based on this observation, the top resizing location can be obtained by finding the first local maximum of intensity summation in each row by the equation S(y)=∑x=1x=XI(x,y) above a threshold (in this case, mean intensity summation is used) to reduce computational cost. Therefore, all rows 15 pixels above are empirically cropped. Similarly, the lower resizing position can be estimated from typical corneal thickness and the axial resolution of the image. The left and right parts are also cropped to remove unwanted structures and the low SNR region on both sides.

As shown in Figure 2, two main types of artifacts are present in OCT images: horizontal artifacts and central noise artifacts. The former appears as long horizontal stripes of high intensity, while the latter is characterized as a vertical saturation region. It is essential to eliminate these artifacts as they have a significant negative impact on the segmentation algorithm.

The horizontal artifact is characterized as adjacent rows with higher mean intensity than others, an efficient way to mitigate it is to subtract the pixel value of each row from the mean intensity of that row.

The central noise artifact is detected by finding a sudden increase in the average intensity of the A-scans, as it is characterized by relatively higher intensity pixels vertically. The whole image is divided equally into three regions and the average intensity (μ) of the A-scans in the peripheral region is calculated. Assuming that the central artifact only occurs in the central region, we therefore consider the A-scan in the central region above a certain threshold [43μ (10)] as the region contaminated by the central artifact. Once the artifact is detected, the region within the artifact is not included in the subsequent processes.

It is necessary to suppress noise in OCT images as the additive thermal and electronic noise can degrade the performance of the algorithm. A 5 × 5 Wiener denoiser (15) is used to increase the image SNR. For simplicity, the image after artifact suppression and denoising is still denoted as I(x, y).

In this section we will focus on approximating the position of the anterior surface boundary, as the air-epithelium interface is generally the region with the best quality (high SNR) in OCT images. The main feature used in the search is the bright-to-dark or dark-to-bright transitions in the axial (vertical) direction. Instead of using a gradient with a directional filter to extract the corresponding boundary (10), a novel adapted estimation method is presented.

It is observed that relatively high pixel intensity occurs at adjacent corneal layer boundaries as a result of over-exposure of reflected and scattered light at the edge in the OCT system. Therefore, the position of the anterior interface can be easily estimated based on the prior assumption that the strongest response of an OCT system (the pixel with the highest intensity) in each A-scan mostly occurs near the corneal boundary instead of random noise. In addition, only pixels with local maxima in A-scans are considered as candidate pixels to reduce computational cost and improve accuracy by excluding other pixels.

To estimate the anterior corneal boundary, we defined a “boundary function” to characterize the corneal interface. In essence, the “boundary function” is a mathematical optimisation objective function that aims to determine the position of the boundary in each A-scan. Therefore, the anterior boundary can be obtained by maximizing the “boundary function” in a predefined region according to pixel resolution and layer thickness to find the optimal input argument. The function is defined in each A-scan in the Equation 1.

Two constraints are used in this function to characterize the corneal boundaries (Equation 2). T⁺(x) is used to account for the difference in pixel intensity between two boundary pixels as a result of over-exposure of scattered and reflected light. The second constraint, T⁻(x), aims to test the lowest intensity of all candidate pixels between two boundary pixels, as there is less scattered light within the corneal layer with similar tissue. The anterior surface is coarsely segmented by finding the top pixel in the axial direction, as shown in Figure 3.

Due to noise in low-SNR regions, the stability problem that some candidate pixels are in random noise should be addressed. To deal with random noise in low-SNR regions (usually in the outer parts), a second-order polynomial approximation (5, 7, 10) is used to fit the profile to eliminate the effect of random noise and form a smooth boundary (other layer boundaries without a predefined approximation model can simply use a median filter to address this problem).

In this section, a novel method is proposed to refine the boundary estimated above. In order to precisely refine the layer boundary, the tactic used is that the latter boundary pixels are determined based on the previously determined boundary pixels with a constant decay weight.

The actual boundary is defined as the maximum intensity change (dark to light or light to dark) in the axial direction. Therefore, the actual layer boundary is found by the maximum absolute vertical gradient. Assuming that the SNR in the center of the corneal image is relatively high compared to the outer part, a non-linear adjustment is considered here, which means that actual boundary pixels at the periphery need more actual central boundary pixels to be confirmed, while central boundary pixels need only a small amount of actual interface pixels to be decided. Therefore, the whole image in the center is divided into two parts to find the actual boundary.

First, the magnitude of the image gradient in the axial direction, symbolized as g(x, y), is calculated using the forward difference gradient operator ∂I∂y=I(x,y+1)-I(x,y)2, and then for each of the latter refined boundary pixels from the center to the periphery is determined based on the maximum summation of the absolute gradient of all previously determined boundary pixels with a geometric distribution decay (constant p) according to Equation 3 by iteratively shifting previous refined pixels up and down in a limited search region (in this case with 5 pixels up and down from the approximated air-epithelium layer interface). The pseudocode for the detailed refinement procedure is shown in Algorithm 1.

After the actual air-epithelium is detected with the proposed method, to smooth the curve of the layer interface, the Savitzky-Golay filter (16–18), with the first-order polynomial and 21-frame length (11), is implemented. Figure 4 shows the refined anterior surface.

Our proposed method for approximating other layers is based on the refined air-epithelium layer interface, as they will have similar boundary profiles. Other layer interfaces can be estimated by maximizing the summation of the absolute gradient, according to the Equation 4, in a region S vertically below the air-epithelium profile (based on typical thickness and pixel resolution). Approximations of other layers in different datasets are shown in Figure 4.

Refinement of other layer interfaces based on the initial estimation in Section 2.5 is performed using a similar technique to that used to adjust the air-epithelium interface in Section 2.4. According to Liu et al. and González-Méijome et al. (19, 20), the cornea is relatively thinner at the center than at the periphery, as shown in Figure 4. Therefore, a normalized linear growth geometric distribution factor p (p∈[0, 1]) from the center to the periphery is used based on the locality of the data points, which means that the approximated data points at the periphery have more ability to explore the layer interface instead of considering the predetermined candidate pixels more than those at the center. However, due to the low SNR at the periphery near the layer boundary, the issue of segmentation stability is raised. To address this issue of refinement flexibility, the refinement technique in the Section 2.4 is again used with a low geometric distribution factor for curve smoothing. Figure 4 shows the final segmentation of all visual corneal layers.

Reconstruction of 3D surface maps of the cornea follows the method for segmentation of 2D B-scans described in the previous sections. There are three main steps: preprocessing, estimation and refinement.

Similar processing steps are implemented to crop the 3D volume image to ensure that only the ROI remains and that various artifacts are removed.

To estimate the air-epithelial interface, the “boundary function” approximation described in Section 2.3 is used with a quadratic surface fit to eliminate the effect of random noise. A similar refinement technique is used in Section 2.4. In 3D segmentation, the starting pixel is in the center of the cornea as we assume that high resolution is presented in the center of the image while low SNR appears in the periphery. For non-linear adaptation, additional information from neighboring pixels was introduced. Candidate pixels are considered using geodesic distance transform (here the city block method is used) (21). Each candidate pixel from the center to the periphery (the distance after transformation) is determined based on the maximum summation of the absolute gradient of all previously refined pixels with geometric distribution decay in a limited region, which means that the algorithm considers more the neighboring candidate pixels and less the distant candidate pixels. A 3 × 3 × 3 median filter is used for curve smoothing.

Other layer interfaces can be approximated by maximizing the sum of the absolute gradient of all estimated pixels in a limited region by shifting the refined air-epithelium profile, following Equation 5.

Then the same refinement method mentioned above is used to find the actual boundary pixels. The final segmentation result is shown in Figure 5.

The algorithm is implemented in Matlab and runs without parallel processing on a Win10 64-bit OS PC with Intel Core i5-7500 CPU @ 3.40 GHz and 8.00 GB RAM.

To evaluate the performance of our method and the other (13), the segmentation results of different layer surfaces were compared with the ground truth described in Section 2.1. The mean unsigned surface position error (MSPE) (13) is used as a metric to evaluate performance. In order to examine the intra- and inter-observer agreement, for Dataset1 the annotations of the same observers and between observers were also compared. For Dataset2, due to data availability, only intra-observer variation was assessed.