To provide a tool for quantifying the effects of retinitis pigmentosa (RP) seen on spectral domain optical coherence tomography images, an automated layer segmentation algorithm was developed. This algorithm, based on dual-gradient information and a shortest path search strategy, delineates the inner limiting membrane and three outer retinal boundaries in optical coherence tomography images from RP patients. In addition, an automated inner segment (IS)/outer segment (OS) contour detection method based on the segmentation results is proposed to quantify the locus of points at which the OS thickness goes to zero in a 3D volume scan. The segmentation algorithm and the IS/OS contour were validated with manual segmentation data. The segmentation and IS/OS contour results on repeated measures showed good within-day repeatability, while the results on data acquired on average 22.5 months afterward demonstrated a possible means to follow disease progression. In particular, the automatically generated IS/OS contour provided a possible objective structural marker for RP progression.
©2011 Optical Society of America
Optical coherence tomography (OCT) has emerged as an important noninvasive imaging modality in the field of ophthalmology, especially with the introduction of spectral domain OCT (SD-OCT) [1–4]. Due to its greater imaging speed and resolution, the layers of the outer retina can be identified on SDOCT scans and the effects of outer retinal diseases, such as retinitis pigmentosa (RP), can be quantified.
RP is a progressive retinal disease that affects the receptors resulting in a severe loss of vision. The structural damage of the outer retina seen on OCT scans has been quantified and used for guidance in therapeutic interventions [5–8]. With SD-OCT, the thickness of the outer segments (OS) of the receptors can be measured and regions of thinning can be identified as well [9–11]. Thinning of the OS layer precedes changes in other receptor layers such as the outer nuclear layer . Recently, Hood et al. argued that the IS/OS contour, the locus of points at which one of the outer retinal borders, the IS/OS line disappears and the OS thickness goes to zero, shows promise as a clinical measure of disease progression . Consequently, it is important to develop and validate automated procedures for segmenting the layers of the outer retina and quantifying the thickness of the OS layer.
Manual segmentation has been employed in previous OCT studies with RP patients [9–12]. However, manual segmentation is a time consuming process and may exhibit subjective variations among different segmentation experts. To provide objective measures and to extend the research study to a larger scale, it is highly desirable to develop a segmentation and quantification methodology that is both reliable and automated. Recently, we have developed a fully automated segmentation algorithm based on gradient information in dual scales, and it showed high accuracy and repeatability in the segmentation of three dimensional (3D) OCT volume scans of normal subjects . Several other algorithms in the literature also have demonstrated automated retinal boundary segmentation based on intensity variation, gradient information, or textural features [14–28]. Among these algorithms, some have shown high accuracy of detecting multiple inner and outer retinal boundaries in normal subjects [18,19,24–26]. Yet the morphological changes in RP diseased eyes that leads to the loss of the IS/OS boundary and thinning of the ONL have remained a technical challenge for the automated segmentation of the OCT images.
The purpose here is to describe modifications to our automated segmentation algorithm  to allow delineation of the inner limiting membrane (ILM) and three outer retinal boundaries in RP OCT images and to provide a validation of this algorithm by comparing the results to those obtained with manual segmentation. In addition, an automated method is proposed to quantify the OS thickness loss in a 3D volume scan and to identify the IS/OS contour in SDOCT images from RP patients.
2.1. Segmentation algorithm
The segmentation algorithm framework is primarily based on the previously reported method , where dual-scale gradient information from a customized Canny edge detector and axial intensity gradient map is employed to create a search graph, and the boundary is then identified by a shortest-path search based on a dynamic programming technique. In this manner, nine inner and outer retinal boundaries are detected as previously described . However, it should be noted that the main focus of the present segmentation is the outer retinal boundaries, particularly the IS/OS border, whose structural changes have been of increasing interest in recent RP OCT studies [8–12]. From this viewpoint, the four boundaries within the scope of this study have been chosen as Boundary 1: between the vitreous and ILM, Boundary 2: between IS and OS, Boundary 3: between the OS and the retinal pigment epithelium (RPE), and Boundary 4: between Bruch’s membrane (BM) and the choroid. (See Fig. 1 .)
The detection order of these four boundaries is the ILM, IS/OS, BM/Choroid and lastly the OS/RPE. The details have been previously reported . In short, the ILM and IS/OS boundaries are detected first so that the detection of the remaining boundaries can be restricted to successively smaller search areas, thereby considerably shortening the computation time. For the detection of each boundary, a graph map based only on node cost assignments is constructed. The node costs are primarily based on a linear combination of Canny edge detector output and axial intensity gradient values. Next the shortest path search algorithm finds the minimum cost path utilizing dynamic programming. In addition, the output from the shortest path search for the OS/RPE and BM/Choroid are smoothed using a polynomial curve-fitting based technique.
Three adjustments to the detection of the BM/Choroid, OS/RPE and IS/OS boundaries have been added to the segmentation algorithm in order to accommodate the corresponding structural changes seen in RP OCT images. After the algorithm detects the IS/OS boundary, the BM/Choroid boundary is detected. The choroidal vascular structure is sometimes more visible in RP OCT images compared to OCT images of healthy controls probably due to the reduced total retinal thickness . Therefore, the high contrast choroidal vascular contour can interfere with the BM boundary detection if the search map were to include only the gradient information, as was the case in the original algorithm . As a result, the search map has been enhanced with the global gradient difference between pixels above the boundary and pixels below the boundary within a base window size of 3 pixels by 8 pixels in the Z- and X-axis directions, respectively. A larger kernel in the X-axis is used because the BM/Choroidal boundary has a lower gradient of change than the contours of the choroidal vasculature. The localized edges of the choroidal vascular contour can be suppressed to some degree via the bigger kernel filter. Figure 2 shows segmentation results derived from both the original BM/Choroid search map (Fig. 2b) and the current enhancement method (Fig. 2c).
Another adjustment is for the OS/RPE detection. Because the OCT image of a RP patient often does not have a visible IS/OS border across a portion of the image, the OS/RPE boundary is allowed to overlap with the IS/OS boundary. The OS thickness is then calculated as the difference between the OS/RPE and IS/OS.
Lastly, the initially detected IS/OS boundary is the dark-to-bright edge of the IS/OS line . As the manual IS/OS segmentation was based on slightly different criteria, the IS/OS boundary is adjusted to detect the location of the area with the highest intensity  rather than the dark-to-bright edge. A second path search favoring high intensities in the image between the initially detected IS/OS and OS/RPE boundaries has been added in order to detect the highest intensity location immediately below the dark-to-bright edge. The graph map used for this secondary shortest path search is generated from the image portion between the initially detected IS/OS boundary (located as the dark-to-bright edge) and the OS/RPE boundary, and the node costs correspond to the inverted pixel intensities. The minimum cost path within the graph is then detected using a dynamic programming algorithm. Lastly, the path is converted back to the original image space to become the final IS/OS boundary favoring high intensity. The result is shown in Fig. 3 .
2.2. IS/OS contour detection method
An automated IS/OS contour detection method for a segmented OCT 3D macular volume is proposed. The IS/OS contour is the locus of points at which the OS thickness goes to zero. The detailed IS/OS contour detection method illustrated in Fig. 4 where the detected IS/OS contour is shown as the yellow contour in Fig. 4d. Firstly, an OS thickness map for a 3D volume is calculated from the distance between the IS/OS and OS/RPE boundaries after the segmentation (Fig. 4a). A morphological filter is then applied to the raw thickness map to remove possible noise (Fig. 4b). The morphological filtering step includes an erosion followed by a dilation operation, both using a disk-shaped filter with a radius of 6 pixels. The thickness map is then converted to binary representation using a threshold of 0 (Fig. 4c), and the largest connected component in the binary image is determined using an 8-connected neighborhood criteria. Finally, the exterior boundaries of the largest connected component are traced using Moore-Neighbor’s tracing algorithm modified by Jacob's stopping criteria  (Fig. 4d). The yellow contour in Fig. 4d shows where the OS thickness becomes zero in this example OS thickness map. The OS quantification algorithm was implemented in Matlab (The MathWorks, Natick, MA, version 2010b).
2.3. Experiment and quantitative analysis
Analyses included the accuracy and repeatability of the segmentation, the repeatability of the IS/OS contour, and the detectability of progression and were performed on the same data from seven patients with RP. Six of these patients participated in a previous study with manual segmentation . The seventh patient exhibited no visible IS/OS border anywhere in the 3D OCT volume, and the data set was used solely in the segmentation accuracy and repeatability study. All patients had SD-OCT three-dimensional macular cube scans (Topcon 3D OCT-1000 system from Topcon Corp., Tokyo, Japan) of both eyes on the first visit. The macular cube scan consisted of 128 B-scans by 512 A-scans with an axial resolution of 3.5 μm/pixel and covered a 6mm by 6mm area. The patients returned for follow-up OCT scans between 10 and 31 months later. At least one eye for all the patients had three-dimensional scans on the follow-up visit (Topcon 3D OCT-2000 system from Topcon Corp., Tokyo, Japan). The macular cube scan also consisted of 128 B-scans by 512 A-scans with an axial resolution of 2.6 μm/pixel and covered a 6mm by 6mm area. Repeat scans were obtained at each session. Written informed consent was obtained from all subjects. Procedures followed the tenets of the Declaration of Helsinki, and the protocol was approved by the Committee of the Institutional Board of Research of Columbia University.
The accuracy of the segmentation of the ILM, IS/OS, OS/RPE and BM/Choroid boundaries was validated by comparing the automated algorithm results to the results with a manual segmentation procedure . The data sets included 88 B-scans taken from 8 RP macular cube scans of 7 patients . Eleven individual B-scans from each cube scan were manually segmented by an experienced segmenter with the use of a computer-aided manual segmentation procedure . The selected B-scans included the line scan centered on the fovea (0°) and 10 others located at: ±1°, ±3°, ±5°, ±7°, and ±9° . The automated segmentation algorithm was applied to the same 3D data sets, and the results for the individually selected B-scans were extracted. The unsigned and signed differences in local border position of the ILM, IS/OS, OS/RPE and BM/Choroid boundaries were calculated . Differences were not calculated for locations where the manual segmenter reported that a boundary was not visible. At each valid A-scan location in the images, the local border position difference was calculated as the difference between the border position for the algorithm and the manual segmentation. The signed and unsigned local border differences were then averaged across A-scans and images.
To test the repeatability of the segmentation algorithm, macular 3D scans of 13 eyes (each with two repetitions) from seven RP patients were segmented. The intra-class correlation (ICC), mean of coefficients of variation (CV) and mean of standard deviations (SD) were calculated to test the repeatability of IS/OS, OS/RPE and BM/Choroid boundaries . For these measures of repeatability, the data within each group of 4 adjacent A-lines were averaged to yield a reduced data set of size 128 by 128. An automated foveal centering algorithm was applied to each repetition for each subject. The data sets were cropped to a center 5mm by 5mm square region, resulting in a final data set size of 107 by 107. The ICC, CV, and SD calculations were then performed between the two repetitions for each subject, and the results from the 13 individual eyes were averaged yielding collective means of ICC, CV, and SD.
After testing the accuracy and repeatability of the segmentation algorithm, the 5mm by 5mm OS thickness maps, and the IS/OS contours detected from two repetitions of the first visit, were compared. A total of 11 macular volumes from 6 subjects were tested. The data sets from the seventh patient were excluded as no OS thickness was detected over the entire volume in each data set from both visits. Two measures were evaluated to quantify the repeat reliability: 1) pair-wise pixel-to-pixel Pearson correlation coefficient was calculated for each pair of 5x5mm OS thickness maps, and 2) the area inside the non-zero OS thickness area contour of the OS thickness maps. In addition, the RP progression detectability of the IS/OS contour was examined in each of the eight eyes among the six patients who had follow-up visits. The OS areas within the IS/OS contours were compared between the average of the two repetitions of the first visit and the average of the two repetitions of the second visit.
3.1. Segmentation results
Figure 5 shows both the manual and automated segmentation results for different degrees of OS thinning. The automated segmentation performs well in locating the peripheral IS/OS boundary disappearance seen on either the nasal or temporal side. Even where the IS/OS border is not visible in the image, the algorithm correctly locates the IS/OS boundary overlapped with the OS/RPE boundary. Figure 6 illustrates typical B-scan segmentation results from a 3D OCT volume of a RP patient eye and the corresponding two-dimensional OS thickness map and IS/OS contour. All data sets were processed by a personal computer (CPU: Quad Core 1.73 GHz, RAM: 4 GB); it took about 10.5 seconds in C++ for the full nine layer detection for each 3D volume. The algorithm has been specifically optimized with multiple threads.
The unsigned and signed border position differences between manual and automatic segmentation are summarized in Table 1 . The unsigned border position difference of the ILM, IS/OS, and OS/RPE boundaries for the algorithm versus the manual segmentation was roughly similar to the axial pixel size of 3.5 µm. All the signed border position differences were less than 1 µm.
The grand average ICC, CV and SD repeatability results for the 13 patient eyes are presented in Table 2 . Overall the ICC of each boundary was 0.98, the mean coefficient of variation was less than 1.28%, and the mean standard deviation was less than 2.77µm. Both the accuracy and repeatability results are comparable to what we have reported on normal subjects .
3.2. OS quantification results
The OS thickness maps and auto-detected IS/OS contours from the two repetitions of the first visit were compared. The first two columns of Fig. 7 show the results from the right eyes of six patients. The pair-wise pixel-to-pixel Pearson correlation coefficient calculated for each pair of 5mmx5mm OS thickness maps was [0.85 0.94 0.97 0.98 0.99 0.98 0.98 0.89 0.88 0.95 0.98] for 11 eyes of these 6 patients. The mean Pearson r was 0.95. The calculated area difference within the IS/OS contour between two repetitions of each eye was graphed as a Bland-Altman plot against the mean OS area of these two repetitions (Fig. 8 ). Most of the difference values are well within the interval limits, that is, average difference ± 1.96 standard deviation of the difference. The outlier in the plot is from the right eye of patient 3. This patient had little OS thickness, often just one pixel (3.5µm) in depth, over a wide area (see third row in Fig. 7). It is difficult for an automatic algorithm to detect the exact two-dimensional contour with high repeatability when the OS thickness is very thin over a relatively broad area. It is also noteworthy that a peninsula-shaped contour was detected in the first visit of that patient (see the images in the first and second columns of the third row in Fig. 7), which probably could not be easily detected via individual B-scans. The IS/OS contour detection methodology proposed here integrates the neighboring B-scan information in two ways. First, the segmentation results were smoothed across neighboring B-scans, and secondly, the OS thickness map was then transformed by a morphological filter before the IS/OS contour was tracked. Therefore, the IS/OS contour derived from a 3D volume based on the above method provides more consistent results across B-scans as compared to the contour based upon the disappearance of the IS/OS line detected solely from each individual 2D B-scans.
Furthermore, the thickness maps of the six patients with follow-up visits (the last two columns in Fig. 7) were compared to those for the first visit. The OS area differences between the average of the two repetitions from the first visit and the average of the two repetitions from the follow-up visit for each of the eight eyes from these six patients are shown in Table 3 , along with the corresponding OS areas from each visit. All the OS area differences are negative, indicating that smaller OS areas were found in the second visit for each patient, and patients 2 and 3 have greater decreases, more than 1mm2, in the OS area. These results show that the OS area calculation may have potential to serve as an indicator for RP disease progression, although it would be beneficial to perform a study with a greater number of patients.
A fully automated segmentation and IS/OS contour detection tool for OCT images from RP patients has been presented and evaluated. The segmentation procedure provides a fast and reliable measurement of the ILM, IS/OS, OS/RPE and BM/Choroid boundaries as indicated by the accuracy and repeatability test results. Furthermore, the IS/OS contour detection tool showed good repeatability and holds promise to assist with RP diagnosis and follow-up treatment. In the future, a larger patient population study should be studied to fully evaluate the automated segmentation and IS/OS contour algorithms.
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(5035), 1178–1181 (1991). [CrossRef] [PubMed]
2. A. F. Fercher, C. K. Hitzenberger, G. Kamp, and S. Y. Elzaizt, “Measurement of intraocular distances by backscattering spectral interferometry,” Opt. Commun. 117(1-2), 43–48 (1995). [CrossRef]
3. M. Wojtkowski, T. Bajraszewski, P. Targowski, and A. Kowalczyk, “Real-time in vivo imaging by high-speed spectral optical coherence tomography,” Opt. Lett. 28(19), 1745–1747 (2003). [CrossRef] [PubMed]
5. M. A. Apushkin, G. A. Fishman, K. R. Alexander, and M. Shahidi, “Retinal thickness and visual thresholds measured in patients with retinitis pigmentosa,” Retina 27(3), 349–357 (2007). [CrossRef] [PubMed]
6. T. S. Aleman, A. V. Cideciyan, A. Sumaroka, E. A. Windsor, W. Herrera, D. A. White, S. Kaushal, A. Naidu, A. J. Roman, S. B. Schwartz, E. M. Stone, and S. G. Jacobson, “Retinal laminar architecture in human retinitis pigmentosa caused by Rhodopsin gene mutations,” Invest. Ophthalmol. Vis. Sci. 49(4), 1580–1590 (2008). [CrossRef] [PubMed]
7. S. G. Jacobson, T. S. Aleman, A. Sumaroka, A. V. Cideciyan, A. J. Roman, E. A. Windsor, S. B. Schwartz, H. L. Rehm, and W. J. Kimberling, “Disease boundaries in the retina of patients with Usher syndrome caused by MYO7A gene mutations,” Invest. Ophthalmol. Vis. Sci. 50(4), 1886–1894 (2009). [CrossRef] [PubMed]
8. S. G. Jacobson, A. V. Cideciyan, T. S. Aleman, A. Sumaroka, E. A. Windsor, S. B. Schwartz, E. Heon, and E. M. Stone, “Photoreceptor layer topography in children with leber congenital amaurosis caused by RPE65 mutations,” Invest. Ophthalmol. Vis. Sci. 49(10), 4573–4577 (2008). [CrossRef] [PubMed]
9. D. C. Hood, C. E. Lin, M. A. Lazow, K. G. Locke, X. Zhang, and D. G. Birch, “Thickness of receptor and post-receptor retinal layers in patients with retinitis pigmentosa measured with frequency-domain optical coherence tomography,” Invest. Ophthalmol. Vis. Sci. 50(5), 2328–2336 (2009). [CrossRef] [PubMed]
10. N. V. Rangaswamy, H. M. Patel, K. G. Locke, D. C. Hood, and D. G. Birch, “A comparison of visual field sensitivity to photoreceptor thickness in retinitis pigmentosa,” Invest. Ophthalmol. Vis. Sci. 51(8), 4213–4219 (2010). [CrossRef] [PubMed]
11. D. C. Hood, M. A. Lazow, K. G. Locke, V. C. Greenstein, and D. G. Birch, “The transition zone between healthy and diseased retina in patients with retinitis pigmentosa,” Invest. Ophthalmol. Vis. Sci. 52(1), 101–108 (2011). [CrossRef] [PubMed]
12. D. C. Hood, R. Ramachandran, K. Holopigian, M. A. Lazow, D. G. Birch, and V. C. Greenstein, “Method for deriving visual field boundaries from OCT scans of patients with retinitis pigmentosa,” Biomed. Opt. Express 2(5), 1106–1114 (2011). [CrossRef] [PubMed]
13. Q. Yang, C. A. Reisman, Z. Wang, Y. Fukuma, M. Hangai, N. Yoshimura, A. Tomidokoro, M. Araie, A. S. Raza, D. C. Hood, and K. Chan, “Automated layer segmentation of macular OCT images using dual-scale gradient information,” Opt. Express 18(20), 21293–21307 (2010). [CrossRef] [PubMed]
14. D. Koozekanani, K. Boyer, and C. Roberts, “Retinal thickness measurements from optical coherence tomography using a Markov boundary model,” IEEE Trans. Med. Imaging 20(9), 900–916 (2001). [CrossRef] [PubMed]
15. 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(6), 2012–2017 (2005). [CrossRef] [PubMed]
16. D. Cabrera Fernández, H. M. Salinas, and C. A. Puliafito, “Automated detection of retinal layer structures on optical coherence tomography images,” Opt. Express 13(25), 10200–10216 (2005). [CrossRef] [PubMed]
17. M. Mujat, R. Chan, B. Cense, B. Park, C. Joo, T. Akkin, T. Chen, and J. de Boer, “Retinal nerve fiber layer thickness map determined from optical coherence tomography images,” Opt. Express 13(23), 9480–9491 (2005). [CrossRef] [PubMed]
18. M. K. Garvin, M. D. Abramoff, R. Kardon, S. R. 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(10), 1495–1505 (2008). [CrossRef] [PubMed]
19. 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(9), 1436–1447 (2009). [CrossRef] [PubMed]
22. V. Kajić, B. Považay, B. Hermann, B. Hofer, D. Marshall, P. L. Rosin, and W. Drexler, “Robust segmentation of intraretinal layers in the normal human fovea using a novel statistical model based on texture and shape analysis,” Opt. Express 18(14), 14730–14744 (2010). [CrossRef] [PubMed]
23. S. Lu, C. Y.-I Cheung, J. Liu, J. H. Lim, C. K.-s. Leung, and T. Y. Wong, “Automated layer segmentation of optical coherence tomography images,” IEEE Trans. Biomed. Eng. 57(10), 2605–2608 (2010). [CrossRef] [PubMed]
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(6), 1321–1330 (2010). [CrossRef] [PubMed]
25. 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(18), 19413–19428 (2010). [CrossRef] [PubMed]
26. V. Kajić, B. Považay, B. Hermann, B. Hofer, D. Marshall, P. L. Rosin, and W. Drexler, “Robust segmentation of intraretinal layers in the normal human fovea using a novel statistical model based on texture and shape analysis,” Opt. Express 18(14), 14730–14744 (2010). [CrossRef] [PubMed]
27. H. Zhu, D. P. Crabb, P. G. Schlottmann, T. Ho, and D. F. Garway-Heath, “FloatingCanvas: quantification of 3D retinal structures from spectral-domain optical coherence tomography,” Opt. Express 18(24), 24595–24610 (2010). [CrossRef] [PubMed]
28. A. Yazdanpanah, G. Hamarneh, B. R. Smith, and M. V. Sarunic, “Segmentation of intra-retinal layers from optical coherence tomography images using an active contour approach,” IEEE Trans. Med. Imaging 30(2), 484–496 (2011). [CrossRef] [PubMed]
29. R. C. Gonzalez, R. E. Woods, and S. L. Eddins, Digital Image Processing Using MATLAB (Pearson Prentice Hall, 2004).
30. D. C. Hood, J. Cho, A. S. Raza, B. A. Dale, and W. Min, “Reliability of a computer-aided manual procedure for segmenting optical coherence tomography scans,” Optom. Vis. Sci. 88, 113–123 (2010). [CrossRef]