Fast retinal layer segmentation of spectral domain optical coherence tomography images

Abstract. An approach to segment macular layer thicknesses from spectral domain optical coherence tomography has been proposed. The main contribution is to decrease computational costs while maintaining high accuracy via exploring Kalman filtering, customized active contour, and curve smoothing. Validation on 21 normal volumes shows that 8 layer boundaries could be segmented within 5.8 s with an average layer boundary error <2.35  μm. It has been compared with state-of-the-art methods for both normal and age-related macular degeneration cases to yield similar or significantly better accuracy and is 37 times faster. The proposed method could be a potential tool to clinically quantify the retinal layer boundaries.


Introduction
Optical coherence tomography (OCT) 1 is an optical signal acquisition and processing modality which is noninvasive for imaging subsurface tissue structures. OCT can attain images of the internal structures of the eye with higher spatial resolution (several micrometers) than other imaging modalities, such as ultrasound, x-ray, and magnetic resonance imaging.
Because human eyes are optically penetrable, OCT has recently become an important technique to revolutionize the clinical imaging of the retina. 2 Quantities derived from the segmented retinal layers will help doctors to understand the anatomy of the human retina better and to quantify the extents of abnormalities. These quantities include nerve fiber layer (NFL), ganglion cell layer (GCL), inner plexiform layer (IPL), inner nuclear layer (INL), outer plexiform layer (OPL), outer nuclear layer (ONL), inner segment (IS), outer segment (OS), and retinal pigment epithelium (RPE) (Fig. 1).
There are efforts to segment layer boundaries from OCT retinal images based on grayscale variation and grayscale thresholding. [3][4][5][6][7][8][9][10] Unfortunately, these methods are sensitive to noise. A Markov random field model for extracting the inner and outer retinal borders was introduced in Ref. 11. Its extension was adopted for optic nerve head 12 and outer retinal layers' segmentation. 13 The autoregressive model was shown to be more robust than those based on grayscale thresholding and grayscale variation. Its main problem was to find reliable "seed" points for OCT images with abnormal retina. In addition, special rules needed to be applied for correcting errors since the model depended on the connection of one-dimensional (1-D) points.
Active contour models were employed to segment retinal layers including macula and optic nerve regions. [14][15][16][17][18] However, they required good initialization and had high computational costs. Methods based on graph theory could accurately segment retinal layer boundaries in normal adult eyes, [19][20][21][22][23][24][25][26][27][28][29] but had high computational complexity. Recently, machine learning approaches, including support vector machines, 30 random forest, 31 and Bayesian artificial neural networks, 32 have attracted much attention. As pointed out in Ref. 10, most of the reported segmentation methods on two-dimensional (2-D) and three-dimensional (3-D) OCT data were not practical for general clinical use due to their high computational costs. However, segmentation will be particularly valuable in fields such as computer-assisted surgery, where real-time visualization of the anatomy is a crucial component. As quantitative changes in different retinal layers are correlated well with changes in visual function and may represent surrogate indicators of glaucoma, 33,34 age-related macular degeneration (AMD), 35 type 1 diabetes, 36 multiple sclerosis, [37][38][39] Alzheimer's disease, 40 and Parkinson's disease, 41 it is of great value to quantify the thickness of different retinal layer boundaries with speed and accuracy. Our objective is to decrease computational costs while preserving reliable and accurate segmentation for spectral-domain OCT (SD-OCT) images. The main points of this study are listed below.
First, Kalman filtering 42 to model the correlation between adjacent image frames is adopted, which, to the best of our knowledge, has not yet been addressed in existing methods. It has the advantages of avoiding the initialization of contours from the second frame and enhancing the robustness in the presence of retinal blood vessel shades and other artifacts of motion or uneven illumination. Second, segmentation of 2-D OCT images is converted to 1-D positioning by exploiting the layered structures of the retina. Finally, a customized active contour model with only an image energy term inspired by Jacob et al. 43 is proposed to roughly localize the retina followed by smoothing based on the Savitzky-Golay algorithms 44 to yield fast and accurate layer boundaries.
The rest of the paper is organized as follows. In Sec. 2, the proposed macular layer boundary segmentation and measurement techniques are presented. Experimental results are then described in Sec. 3, whereas discussion and conclusion are given in Sec. 4.

Proposed Method
An OCT imaging system can obtain 2-D images (i.e., B-scans, where the number of A-scans is the width of a 2-D image) and 3-D volumes (volumetric data, containing spatially aligned multiple B-scans). Due to the information correlation among adjacent images, for volumetric data, the first frame and others are addressed separately. If a volumetric data is not available, the algorithm to process individual B-scans is the same as that to process the first frame of a volumetric data. We will divide Sec. 2 into eight subsections, which, respectively, deal with image projection to get the region of interest (ROI), directional filtering, multiresolution to estimate vitreous-NFL and IS/OS, customized active contour model for coarse segmentation, curve smoothing, image flattening and flatttening, Kalman filtering, and segmentation of other layer boundaries. The following flowcharts outline the main steps (Fig. 2), which are, respectively, for individual B-scans or the first frame of a volumetric data, and the nonfirst frame of a volumetric data.

Image Projection to Get Region of Interest
In OCT images, the data are substantially redundant in the axial direction. Searching for the ROI is essential to decrease the computational cost as well as to make the delineation more reliable. Once the ROI is determined, the subsequent processings will be performed within the ROI. A method of image projection is employed to attain the ROI.
First, each row of the B-scan to be processed is projected along the lateral direction: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 3 2 6 ; 4 2 8 where fðx i ; y i Þ is the grayscale of pixel ðx i ; y j Þ, i and j are, respectively, the horizontal (lateral direction) position and vertical (axial direction) position, and M is the width of the image. The highest and the second highest peaks correspond, respectively, to the center line of the RPE-IS/OS complex and the  center line of the NFL. Here, the center line is a horizontal line whose Y coordinate is the average Y coordinates of the corresponding layer boundaries. Second, according to prior knowledge, the distance between the center line of RPE-IS/OS complex (center line of NFL) and the top (bottom) border of the ROI can be determined, which has a correlation with the height of the B-scan. That is, the distance can be set as ρH. Here, ρ is a predefined constant and H is the height of the B-scan. In this study, we set ρ ¼ 0.1. Figure 3 shows the derived ROI (the rectangle formed by the magenta and yellow lines).
Note that image projection to derive the ROI is only applied for processing individual B-scans or the first frame of a volumetric data. For the nonfirst frame of a volumetric data, Kalman filtering (see Sec. 2.7) is utilized to avoid redundant calculation of image projections frame by frame.

Directional Gaussian Filtering
The original B-scan is smoothed with 1-D Gaussian filter along the lateral direction as the layer boundaries are almost horizontal to suppress noise without blurring layer boundaries. Figure 4 shows the original and smoothed images along the lateral direction with the standard deviation (SD) σ being 5.0 and the mean being 0.

Multiresolution to Estimate Vitreous-Nerve Fiber
Layer and Inner Segment/Outer Segment The two most outstanding layer boundaries are the vitreous-NFL and IS-OS in an SD-OCT retinal B-scan. As such, these two boundaries could be extracted first. To be robust to noise and reduce computational costs, a multiresolution approach is adopted in four steps. First, downsampling is performed. Suppose that the original image is denoted as f H 0 W 0 ðx; yÞ, then the downsampling is carried out by a factor of 2 vertically to get f H 1 W 0 ðx; yÞ. E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 3 2 6 ; 6 1 7 This downsampling could be performed along both the X and Y directions: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 3 2 6 ; 5 1 2 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 3 2 6 ; 4 7 4 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 6 3 ; 5 0 2 f 42 ðx; yÞ ¼ f H 4 W 2 ðx; yÞ Second, the grayscale gradient magnitude gðx; yÞ is calculated at the coarse scale.
E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 6 3 ; 3 0 8 gðx; yÞ ¼ ½g 2 x ðx; yÞ þ g 2 y ðx; yÞ 1∕2 : Third, the two initial contours are composed of pixels whose gradient magnitudes are maximum and the second maximum at each column.
Finally, the edge pixels at the coarse scale are converted to pixels of the original image by multiplying x coordinates by 4 and y coordinates by 16 to form the initial contours (Fig. 5).
Note that multiresolution for estimating vitreous-NFL and IS/OS is only employed for individual B-scans or the first frame of a volumetric data. For the nonfirst frame of a volumetric data, Kalman filtering (see Sec. 2.7) is employed instead.

Customized Active Contour Model for Coarse Positioning of Layer Boundaries
A customized active contour model is proposed to cater to the layered structures in OCT macular images. An active contour [45][46][47] is traditionally modeled as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 3 2 6 ; 3 8 4 where αðsÞ, βðsÞ, and γðsÞ are weights. E continuity is the continuity term to encourage even spacing of the contour points and could be calculated by v i − v i−1 with v i for the i'th contour point. E curvature is the curvature term to penalize rapid curve bending and could be calculated by jv i−1 − 2v i þ v iþ1 j. E image is the image term for the contour to be converged to edges and could be calculated as ðg cur − g min Þ∕ðg max − g min Þ, with g cur , g max , and g min being, respectively, the gradient magnitude of the current pixel, and the maximum and minimum  gradient magnitudes of the current neighborhood. In this study, the following constraints are imposed for each layer boundary: 1. There is only one contour point at each column.
2. At each iteration, the contour point at each column could only move vertically, being either one pixel up or down.
3. The contour point at each column should be adjacent to the contour point of the previous and next columns (eight neighbors).
With these constraints, there might be only four configuration types for the three neighboring contour pixels (Fig. 6).
Here are some observations to justify the simplification of Eq. (11).
1. There are only two cases between the two adjacent contour points with E continuity being either 1 or 2 1∕2 .
2. The curvature term E curvature is 0 for the horizontal and diagonal cases, 1 for the inflectional case, and 2 for the V-shaped case.
Equation (11) could thus be simplified to include only the third term. In this study, the third term is defined locally to reduce the computational cost.
where k is the position of the current pixel and n is the neighborhood size in the vertical direction. Figure 7 shows the derived initial contours of vitreous-NFL and IS/OS from Eq. (12) of Fig. 5.

Curve Smoothing
As the customized active contour model [Eq. (12)] has not taken into account continuity and curvature terms, the derived contour could be jagged. To make it smooth, curve smoothing is employed. Savitzky-Golay smoothing filters 44 have been widely used for 1-D filtering. The fundamental idea is to fit a polynomial to the data surrounding each data point. The difference between the general least-squares polynomial smoothing filters and the Savitzky-Golay smoothing filters lies in the fact that the latter introduces tables of convolution weights for the center-point, which could substantially enhance the speed. 48 In this study, Savitzky-Golay smoothing filters are employed for curve smoothing of the layer boundaries with the order of the polynomial being 4 and a neighborhood size from 20 to 30 pixels for different layer boundaries.

Image Flattening and Unflattening
To handle the possible deformation of layer boundaries due to pathology, the image is flattened with respect to the extracted layer boundary of IS/OS. Specifically, the flattening process is carried out this way: find the maximum y coordinate of the IS/OS layer boundary (denoted as max y ), for each column j, the column will move (max y − y j ) downwards with y j being the y coordinate of the IS/OS boundary pixel. In this way, the extracted IS/OS becomes a horizontal line. Figure 8 shows the flattened image from Fig. 9. The other six layer boundaries are all extracted from the flattened image. The 8 layer boundaries are then derived by reversing the flattening process. It should be noted that these flattening and unflattening procedures are only applied to individual B-scans or the first frame of a volumetric data.

Kalman Filtering
As a prior knowledge, for a volumetric data, the positions of layer boundaries in adjacent images (frames) are similar. The coarse layer boundaries of the current frame could be estimated from the layer boundaries of the previous frame. We employ Kalman filtering 42 to track coarse layer boundaries frame by frame.
To approximate the layer boundaries, Kalman filter is formulated as a state equation [Eq. (13)] and a measurement equation [Eq. (14)]: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 6 3 ; 2 8 9 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 6 3 ; 2 4 7 z k ¼ Hϕ k þ v k ; where ϕ ¼ y v y is a 2-D vector with y being the Y coordinate of layer boundaries and v y being the speed of adjacent frames of the layer boundaries being tracked, w y the process noise, v k the measurement noise, k the current frame, and k − 1 the previous with dt being the time interval from the previous frame to the current frame. B is the control input matrix and will not take effect when there is no control input (here, we suppose u ¼ 0). As there is no control input, the velocity cannot be measured, so z ¼ z y , and The state ϕ and variance P could be predicted as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 3 2 6 ; 3 9 4 where − is for the predicted value and T is for transpose. The correction equation is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 3 2 6 ; 3 4 6 where R is the SD of the measurement noise. K will be adaptive to R, i.e., the components of K will be smaller when R is large, and larger otherwise. Then, the state ϕ can be updated E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 8 ; 3 2 6 ; 2 7 0 At last, the variance P can be updated E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 9 ; 3 2 6 ; 2 2 7 To summarize, two phases [prediction by Eqs. (15) and (16) and updating by Eqs. (17)- (19)] of Kalman filtering are applied to attain the state with minimum variance.

Other Boundaries Segmentation
Once the two most outstanding layer boundaries are segmented, they can be utilized as constraints to facilitate segmentation of other layer boundaries. Figure 10 illustrates the flowchart of segmentation for other layer boundaries.
The segmentation of other layer boundaries consists of the same three steps: initialization, coarse segmentation by the customized active contour method [Eq. (12)], and curve smoothing through the Savitzky-Golay filtering described in Sec. 2.5.  The second way of initialization is for the nonfirst frame of a volumetric data using the Kalman filter described in Sec. 2.7, by employing the fact that layer boundaries between consecutive frames are similar.

Experiments
In this section, experiments were carried out to validate the proposed method. Validation was first carried out on 21 SD-OCT volumes with each volume from one subject was collected from 21 healthy Chinese subjects and then comparisons were made with state-of-the-art methods for both normal and pathological eyes with AMD using common datasets. Finally, additional experiments were performed to vary the components of the algorithm for comparison.
For each layer of a B-scan, the thickness between neighboring boundaries was computed and then the average difference in layer thickness between the automatic and manual delineation was calculated for quantification. The average difference in layer thickness between two manual delineations was also calculated to reflect the possible uncertainty of measuring the layer thickness, even when done by experts. Figure 11 showed the segmentation for a B-scan of a normal right eye. Figure 12 illustrated that the proposed method could handle images with heavy noise and blood vessel artifacts.

Validation on Our Clinical Data
A customized SD-OCT system was adopted for 3-D imaging, which used a superluminescent diode with a full-width at half-maximum spectral bandwidth of 50 nm and a center wavelength of 840 nm, and a linear CCD camera (Aviiva-SM2-CL-2014, e2V) with a 14-μm pixel size operating in an 11-bit mode.   The optical power of the beam on the cornea was not greater than 750 μW. Images were obtained with a depth of 1.68 mm and axial resolution of 3.5 μm.
A dataset consisting of 21 SD-OCT volumes (128 × 2048 × 2048 voxels, 6 × 6 × 2.5 mm 3 macular regions, and an in-plane pixel size of 3.5 μm) with each volume from one subject was collected from 21 healthy Chinese subjects at the Xinhua Hospital Affiliated to Shanghai Jiaotong University School of Medicine. One eye from each subject was selected randomly. The protocol of the research has been approved by the Institutional Review Board of the hospital. All subjects gave written consent and provided permission for educational and scientific purposes.
The automated detection of macular layer boundaries was performed successfully on all the 21 volumes. For segmenting 8 layer boundaries, it took an average of 5.8 s for each volume (128 B-scans per volume) and 95 ms for a B-scan on a personal computer (C program running on a 3.2 GHz Intel core i5-3470 CPU). These data were also segmented manually by two eye specialists as the reference for comparison with automatic segmentation.
Here ×8 was used to process eight boundaries for each B-scan. For processing a volumetric data without the first frame, the time (5715 ms) distribution was: 1-D Gaussian filter (5 ms × 127 ¼ 635 ms), Kalman filtering (1 ms × 8 × 127 ¼ 1016 ms), coarse segmentation (2 ms × 8 × 127 ¼ 2032 ms), and curve smoothing (2 ms × 8 × 127 ¼ 2032 ms). Here ×127 was used to process 127 B-scans of each volumetric data without the first frame. Table 1 showed the absolute mean difference and SD of the 21 volumes.

Comparison with a State-of-the-Art Method for Normal Eyes
This is to compare the proposed method with the method of Chiu et al. 23 using their data and their manual segmentation as a reference for normal eyes. The 10 volumes were acquired by Bioptigen, Inc. (Research Triangle Park, North Carolina) imaging systems with an axial spacing of 3.23 μm. Only 100 B-scans from the 10 volumes were manually delineated by two experts (one junior and one senior) and the manual delineation by the senior expert was taken as the reference, among which 29 were used for interexpert comparison.
The average times to derive the 8 layer boundaries of a B-scan for the proposed method and Chiu's method 23 were, respectively, 60 and 9740 ms.
The comparison results between two experts as well as between the proposed method with Ref. 23 were summarized in Table 2. Paired t-tests were carried out to find that there existed a significant difference in the accuracy of INL (p < 0.004), OPL (p < 0.007), ONL-IS (p < 0.019), OS (p < 0.019), RPE (p < 0.0001), and total retina (p < 0.019), while errors of the NFL and GCL-IPL did not have significant differences.

Comparison with a State-of-the-Art Method for Eyes with Age-Related Macular Degeneration
AMD is a primary cause of vision loss worldwide. 49 To validate the applicability of the proposed method to data with pathology, the second direct comparison was to compare the proposed method with the method of Chiu et al. 26 using their data and their manual segmentation as a reference for AMD patients.  In the study of Chiu et al. 26 for SD-OCT images of eyes containing drusen and geographic atrophy (GA), segmentation was only performed for three retinal layer boundaries [the inner internal limiting membrane (ILM), inner RPE drusen complex, and outer Bruch's membrane defined in Ref. 26 corresponded to, respectively, vitreous-NFL, OS-RPE, and RPE-choroid in Fig. 1]. Figure 13 showed the segmentation of a B-scan of an AMD patient from Ref. 26. In the Chiu's study, 220 B-scans across 20 patients (one volume per patient) for eyes with AMD were acquired by Bioptigen, Inc. (Research Triangle Park, North Carolina) imaging systems with an axial pixel resolution ranging from 3.06 to 3.24 μm per pixel (Devers Eye Institute, 3.21 μm; Duke Eye Center, 3.23 μm; Emory Eye Center, 3.06 μm; and the National Eye Institute, 3.24 μm). All the SD-OCT images were manually delineated by two certified experts (one junior and one senior) and the manual delineation by the senior expert was taken as the referencel all of the images were used for interexpert comparison.
The comparison results between the proposed and Ref. 26 were summarized in Table 3. Paired t-tests were carried out to find that there existed a significant difference in the accuracy of: total retina of drusens with low quality (p < 0.004), total retina of GA with high quality (p < 0.0001), total retina of GA with low quality (p < 0.032), and total retina for both the drusens and GA images (p < 0.0001), while errors of the rest did not have significant differences.

Additional Experiments
To see the comparative performance with respect to different initializations for a volumetric data from the second frame (with the subsequent processing steps unchanged), two variants were devised: one variant to use the segmented layer boundaries of the previous frame as the corresponding initial layer boundaries of the current frame and this variant was denoted as method A; and the other variant was to ignore the spatial correlation to derive initial layer boundaries independently in the same way as was done for the first frame, and this variant was denoted as method B. Table 4 summarized the performance of the original method and the variants.
In refining the boundary layers from coarse segmentation, the proposed method employed a combination of the customized/simplified active contour model and Savitzky-Golay smoothing filter. A variant was to use the full-active contour model. To this end, the active contour model in Ref. 46 was implemented for comparison, and this variant was denoted as method C. The comparative performance was summarized in Table 5.

Sensitivity to Parameters
The parameter ρ determines the size of the ROI. We alter the ρ in the range from 0.05 to 0.40 to yield similar segmentation accuracy (with the total retinal thickness showing a difference of 0.43 AE 0.52 μm), but different time consumption (the average time is in the range from 71 to 116 ms). A ρ value of 0.1 seems a good trade-off between segmentation accuracy and time consumption.
The parameter σ is related to the smoothing degree of Gaussian filtering. As changing the σ in [1.0, 10.0] yields almost the same segmentation, a value of 5.0 is chosen for this study.

Comparison with Existing Methods
The performance of existing methods is to be compared in terms of computational cost and accuracy.
Segmentation of nerve fiber for a single image (with 1000 A-scans) could take as long as 62 s. 14 It took 24 s to segment seven layer boundaries for a 1024 × 512 pixels' image. 5 The RPE segmentation took 8.3 min for a volume with 60 images (1024 × 1000 pixels). 50 As to processing normal OCT volumetric data, the segmentation time for RPE and vitreous-NFL on 320 × 1024 × 138 volumetric data was 37 s, 10 and for 11 layer boundaries on 200 × 1024 × 200 volumetric data was 70 s. 22 The average computation time was 107 s∕volume (400 × 1000 × 11 voxels) for 8 layer boundaries. 23 The fastest approach was 16 s for detecting nine layer boundaries per volume (480 × 512 × 128 voxels) by a commercial system named Topcon 3-D OCT-1000 (Topcon Medical Systems, Oakland, New Jersey). 24 The proposed method segments 8 layer boundaries within 5.8 s for a volumetric data and 95 ms for an individual B-scan with similar spatial resolutions (128 × 2048 × 2048 voxels), being substantially faster than these existing methods. Garvin et al. 19,20 reported an overall mean absolute boundary differences of 6.10 AE 2.90 μm 19 and 5.69 AE 2.41 μm. 20 Better accuracies were achieved to be in the range of 3 to 4 μm, i.e., 3.39 AE 0.96 μm, 24 3.04 AE 2.65 μm, 23 and 3.38 AE 4.09 μm. 31 The proposed method yields a comparable maximum absolute boundary difference of 2.35 AE 1.32 μm ( Table 1). As these accuracies are from the segmentation of different data, it is inherently difficult to judge which method yields the best accuracy, as the segmentation accuracy depends on other factors in addition to segmentation method, such as imaging type [timedomain OCT (TD-OCT), or SD-OCT], the imaging quality, and axial and lateral resolutions.
We have tried to make direct comparisons using common datasets with existing methods.
On the common dataset with normal eyes, 23 the proposed method has been compared with Chiu's method. 23 From Table 2, it can be seen that the proposed method yields similar accuracy to that of Ref. 23  On the common dataset of patients with AMD, 26 the proposed method has been compared with Chiu's method. 26 From Table 3, it could be seen that the proposed method yields similar segmentation accuracy to that of Ref. 26 for the pathological data (9 out of 10 quantities of the proposed method have a smaller average error than Chiu's method 26 ). For the segmentation of three layer boundaries of 110 B-scans with drusens and 110 B-scans with GA, 26 the proposed method has a significantly higher accuracy than Ref. 26 for the total retina of drusens with low quality, total retina of GA with high quality, total retina of GA with low quality, and total retina for both the drusens and GA images, while errors of the rest do not have significant differences. For these pathological data, the proposed method is 37 (1700/45) times faster than Chiu's method 26 to derive the three layer boundaries of a B-scan.

Advantages
The recent evolution of OCT from the time domain to the spectral domain has substantially increased imaging speed and made the acquisition of volumetric data with a huge number of bytes feasible. Therefore, fast and automatic segmentation of layers and thickness measurement for OCT data has become increasingly more important than ever. Our purpose is to minimize processing time while maintaining reliable segmentation results. The proposed method has the following advantages.
1. Fast. The average time cost of the proposed method to extract 8 layer boundaries is 5.8 s for volumetric data and 95 ms for a B-scan, while existing methods will take several minutes. 10 The low computational cost can be attributed mainly to the framework of the proposed method: image projection to initialize the two most prominent layers (RPE-IS/OS and NFL), 1-D Gaussian filtering along lateral direction instead of 2-D Gaussian filtering, approximating the RPE-IS/ OS and NFL based on grayscale gradient magnitude at low spatial resolution (one 16th of the original Y and quarter of the original X), customized active contour model (only calculating grayscale difference along the lateral direction) to refine the layer boundaries, Savitzsky-Golay smoothing with precalculated weights, and Kalman filtering with a low order (2 × 2) matrix manipulation to initialize layer boundaries from previous frame for volumetric data.
2. Accurate. For both normal eyes and eyes with AMD, the proposed method yields similar or higher accuracy as compared with the state-of-the-art methods 23 ( Table 2) and Ref. 26 (Table 3). This could be ascribed to the core components of the proposed method: employment of prior knowledge to derive ROI, Kalman filtering to perform well even in the presence of heavy noise and artifacts, and customized active contour plus Savitzsky-Golay smoothing to balance between converging to edges and smoothness.
3. Being able to handle pathologies. According to our methodology, the pathology will only influence the initialization of layer boundaries of the first frame of a volumetric data, as the initialization of layer boundaries of subsequent fames will be cast from the previous frame. Due to their edge prominence, the virtreous-NFL and IS/OS could be extracted irrespective of pathologies. To get rid of the influence of a possible deformation of layer boundaries due to pathologies for any individual B-scan or the first frame of a volumetric data, the image is flattened such that the extracted IS/OS boundary is a horizontal line, and the other six layer boundaries will be segmented from the flattened image by taking the horizontal line of IS/OS as the first basis with subsequently extracted layer boundaries being possibly used as bases for those layer boundaries to be extracted from. The method has been validated on pathological eyes with drusens and GA to perform well ( Fig. 13 and Table 3).

Contributions and Limitations
From Table 4, it can be found that initialization of layer boundaries through Kalman filtering could yield more accurate segmentation than the direct casting and independent initialization. As a matter of fact, this feature of Kalman filtering may be employed for initialization of other structures of volumetric data, such as the brain from 3-D computed tomography or magnetic resonance images once the brain from a 2-D slice is available.
From Table 5, it can be seen that the proposed combination of a customized/simplified active contour model and Savitzky-Golay smoothing is better than the active contour method 46 in terms of both segmentation accuracy and computational cost. This feature may be explored as one more option to replace the active contour methods.
The novelties of the study are fourfold: determination of the ROI from grayscale projection and prior knowledge, Kalman filtering to initialize layer boundaries from the previous frame, customized active contour for coarse delineation of layer boundaries, and curve smoothing with precalculated weights. All these contribute to substantial reduction of computational cost while preserving segmentation accuracy.
The study is not without limitations. In particular, this study is focused on methodology and validation has only been carried out for retinal layer boundaries in OCT images of normal eyes and eyes with AMD. In the future, we will be handling more pathological cases.
The proposed method is based on the assumption that the layer boundaries are visible in the images and different layer boundaries will not intersect. We humbly believe that the proposed method will be applicable to those pathologies as long as they will not invalidate this assumption.
Existing methods could segment different numbers of layer boundaries (and/or surfaces), i.e., 2 boundaries (Ref. 9 with SD-OCT for posterior retinal layers of human eyes, 10 with SD-OCT for AMD of human eyes, 11 with SD-OCT for total retinal thickness of human eyes, 12 with TD-OCT for total retinal thickness of human eyes, 13 with TD-OCT for optic nerve head of human eyes, 14 with SD-OCT for NFL thickness of human eyes, 15 with SD-OCT for optic nerve head of human eyes), 3 boundaries with SD-OCT for retinal layers of human eyes), varying from OCT types (TD-OCT or SD-OCT), purposes (total retinal thickness, optic nerve head, intraretinal, outer retinal, or pathology), and subjects (human or animal). Among the methods segmenting more layer boundaries than the proposed method, Refs. 17 and 28 focused on rat (or mouse) retina, which are different from human retina; as to Refs. 22, 24, and 31, they adopted Cirrus HD-OCT machines (Carl Zeiss Meditec, Inc., Dublin, California), Topcon 3-D OCT-1000 equipment (Topcon Medical Systems, Oakland, New Jersey), and Spectralis OCT system (Heidelberg Engineering, Heidelberg, Germany) respectively, which are all commercial systems with better signal-tonoise ratio than our customized system. In this study, as experts are only able to delineate 8 layer boundaries due to image quality, we opt to segment these layer boundaries to demonstrate the effectiveness of the methodology. We believe that the proposed methodology could be easily extended for extraction of more layers when the imaging quality is enhanced to be comparable with that of existing commercial systems.
To conclude, in this paper, an approach for automatic measurement of macular thickness using SD-OCT has been proposed and validated. A quantitative evaluation has been performed on 21 volumetric data covering the usually observed variability, by comparing the results with those obtained manually by two experts. For segmenting 8 layer boundaries on each B-scan, quantitative results show that the average layer boundary error is below 2.35 μm, being smaller than the difference between two eye experts' manual delineations with an average processing time of 5.8 s for a volumetric data (128 × 2048 × 2048 voxels) and 95 ms for an individual B scan. The algorithm has been compared with state-of-the-art methods (Refs. 23 and 26) on a common dataset of normal eyes 23 and a common dataset of AMD patients 26 to yield similar or higher accuracy (Tables 2  and 3) and is 37 times faster to derive the 8 or 3 layer boundaries of a B-scan. The proposed method could be a potential tool to quantify the retinal layer boundaries of normal subjects and AMD patients from SD-OCT. In the future, we are going to extend the methodology to more layer boundaries on one hand, and to more eye diseases, such as diabetic macular edema, 51 for computer-assisted diagnosis on the other hand, from SD-OCT images.
Tianqiao Zhang is a PhD student with Shenzhen Institutes of Advanced Technology since 2011. He received his master's degree from Guangxi Normal University of China in 2006. Then he worked as a software engineer for Huawei Technologies Co., Ltd., a senior image analysis engineer, and a team leader with Shenzhen Moptim Imaging Technique Co., Ltd. His research interests include medical image processing and machine vision.
Zhangjun Song received his PhD from Tsinghua University in 2007. Then he worked as an assistant researcher, was promoted to and has been working as a principal investigator with Shenzhen Institutes of Advanced Integration Technology, Chinese Academy of Sciences. He also worked as a research fellow in National University of Singapore for more than one year. His research interests include machine vision and motion control of mobile robots.
Xiaogang Wang has been with the Shanxi Eye Hospital since 2014 as a resident of ophthalmology. He obtained the MD degree from Shanghai Jiao Tong University. His research activities are mainly about ophthalmic imaging technologies for the eye, especially about optical coherence tomography, and Scheimpflug imaging technology.
Huimin Zheng is a research assistant at Shenzhen Institutes of Advanced Technology. She received her BS degree in information and computing technology from Guilin University of Electronic Technology in 2010. Her current research interest is medical image processing.
Fucang Jia received his PhD in computer application technology from Institute of Computing Technology, Chinese Academy of Sciences in 2004. He then worked as an engineer and senior engineer in Shenzhen Anke High-Tech Co., Ltd. From 2008, he joined Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, as one principle investigator. His research interests include medical image processing and computer assisted surgery.
Jianhuang Wu received his PhD from Shenyang Institute of Automation, Chinese Academy of Sciences, in 2007. He is currently an associate professor at the Research Laboratory for Medical Imaging and Digital Surgery, Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences. His research interests include virtual surgery, medical visualization, and computer graphics.
Guanglin Li received his PhD from Zhejiang University in 1997. He was a research associate at University of Illinois at Chicago. Then he joined the BioTechPlex Corporation as a senior scientist, and later became a senior research scientist in the Neural Engineering Center for Artificial Limbs at the Rehabilitation Institute of Chicago. He is currently a professor with Shenzhen Institutes of Advanced Technology. His current research interests include neural engineering and computational biomedical engineering.
Qingmao Hu is a professor with Shenzhen Institutes of Advanced Technology since 2007. He received his PhD from Huazhong University of Science and Technology in 1990. Then he worked as a lecturer and associate professor for the South Medical University of China, a guest scientist with University of Bern in Switzerland, and a senior scientist for A*STAR in Singapore. His research interests include medical image processing, pattern recognition, and computer-aided diagnosis.