Automated segmentation and characterization of choroidal vessels in high-penetration optical coherence tomography

An automated choroidal vessel segmentation and quantification method for high-penetration optical coherence tomography (OCT) was developed for advanced visualization and evaluation of the choroidal vasculature. This method uses scattering OCT volumes for the segmentation of choroidal vessels by using a multi-scale adaptive threshold. The segmented choroidal vessels are then processed by multi-scale morphological analysis to quantify the vessel diameters. The three-dimensional structure and the diameter distribution of the choroidal vasculature were then obtained. The usefulness of the method was then evaluated by analyzing the OCT volumes of normal subjects. © 2013 Optical Society of America OCIS codes: (170.4470) Ophthalmology; (170.4500) Optical coherence tomography; (170.3890) Medical optics instrumentation; (100.2960) Image analysis; (110.4500) Optical coherence tomography. References and links 1. R. A. Linsenmeier and L. Padnick-Silver, “Metabolic dependence of photoreceptors on the choroid in the normal and detached retina,” Invest. Ophthalmol. Vis. Sci. 41, 3117–3123 (2000). 2. D. L. Nickla and J. Wallman, “The multifunctional choroid,” Prog. Retin. Eye Res. 29, 144 – 168 (2010). 3. R. J. Klein, C. Zeiss, E. Y. Chew, J.-Y. Tsai, R. S. Sackler, C. Haynes, A. K. Henning, J. P. SanGiovanni, S. M. Mane, S. T. Mayne, M. B. Bracken, F. L. Ferris, J. Ott, C. Barnstable, and J. Hoh, “Complement factor h polymorphism in age-related macular degeneration,” Science 308, 385–389 (2005). 4. S. Alam, R. J. Zawadzki, S. Choi, C. Gerth, S. S. Park, L. Morse, and J. S. Werner, “Clinical application of rapid serial fourier-domain optical coherence tomography for macular imaging,” Ophthalmology 113, 1425– 1431 (2006). 5. V. Manjunath, J. Goren, J. G. Fujimoto, and J. S. Duker, “Analysis of choroidal thickness in age-related macular degeneration using spectral-domain optical coherence tomography,” Am. J. Ophthalmol. 152, 663–668 (2011). 6. D. R. Guyer, L. A. Yannuzzi, J. S. Slakter, J. A. Sorenson, A. Ho, and D. Orlock, “Digital indocyanine green videoangiography of central serous chorioretinopathy.” Arch. Ophthalmol. 112, 1057–1062 (1994). 7. L. A. Yannuzzi, K. T. Rohrer, L. J. Tindel, R. S. Sobel, M. A. Costanza, W. Shields, and E. Zang, “Fluorescein angiography complication survey.” Ophthalmology 93, 611–617 (1986). 8. M. Hope-Ross, L. A. Yannuzzi, E. S. Gragoudas, D. R. Guyer, J. S. Slakter, J. A. Sorenson, S. Krupsky, D. A. Orlock, and C. A. Puliafito, “Adverse reactions due to indocyanine green.” Ophthalmology 101, 529–533 (1994). 9. A. F. Fercher, W. Drexler, C. K. Hitzenberger, and T. Lasser, “Optical coherence tomography — principles and applications,” Rep. Prog. Phys. 66, 239–303 (2003). #188755 $15.00 USD Received 12 Apr 2013; revised 9 Jun 2013; accepted 10 Jun 2013; published 25 Jun 2013 (C) 2013 OSA 1 July 2013 | Vol. 21, No. 13 | DOI:10.1364/OE.21.015787 | OPTICS EXPRESS 15787 10. W. Drexler and J. G. Fujimoto, “State-of-the-art retinal optical coherence tomography,” Prog. Retin. Eye Res. 27, 45–88 (2008). 11. R. F. Spaide, H. Koizumi, and M. C. Pozonni, “Enhanced depth imaging spectral-domain optical coherence tomography,” Am. J. Ophthalmol. 146, 496 – 500 (2008). 12. R. Margolis and R. F. Spaide, “A pilot study of enhanced depth imaging optical coherence tomography of the choroid in normal eyes,” Am. J. Ophthalmol. 147, 811 – 815 (2009). 13. A. Unterhuber, B. Považay, B. Hermann, H. Sattmann, A. Chavez-Pirson, and W. Drexler, “In vivo retinal optical coherence tomography at 1040 nm enhanced penetration into the choroid,” Opt. Express 13, 3252–3258 (2005). 14. E. C. Lee, J. F. de Boer, M. Mujat, H. Lim, and S. H. Yun, “In vivo optical frequency domain imaging of human retina and choroid,” Opt. Express 14, 4403–4411 (2006). 15. Y. Yasuno, Y. Hong, S. Makita, M. Yamanari, M. Akiba, M. Miura, and T. Yatagai, “In vivo high-contrast imaging of deep posterior eye by 1-um swept source optical coherence tomography andscattering optical coherence angiography,” Opt. Express 15, 6121–6139 (2007). 16. L. Duan, M. Yamanari, and Y. Yasuno, “Automated phase retardation oriented segmentation of chorio-scleral interface by polarization sensitive optical coherence tomography,” Opt. Express 20, 3353–3366 (2012). 17. T. Torzicky, M. Pircher, S. Zotter, M. Bonesi, E. Götzinger, and C. K. Hitzenberger, “Automated measurement of choroidal thickness in the human eye by polarization sensitive optical coherence tomography,” Opt. Express 20, 7564–7574 (2012). 18. R. Motaghiannezam, D. M. Schwartz, and S. E. Fraser, “In vivo human choroidal vascular pattern visualization using high-speed swept-source optical coherence tomography at 1060 nm,” Invest. Ophthalmol. Vis. Sci. 53, 2337–2348 (2012). 19. Y.-J. Hong, S. Makita, F. Jaillon, M. J. Ju, E. J. Min, B. H. Lee, M. Itoh, M. Miura, and Y. Yasuno, “Highpenetration swept source doppler optical coherence angiography by fully numerical phase stabilization,” Opt. Express 20, 2740–2760 (2012). 20. F. Jaillon, S. Makita, and Y. Yasuno, “Variable velocity range imaging of the choroid with dual-beam optical coherence angiography,” Opt. Express 20, 385–396 (2012). 21. B. Braaf, K. A. Vermeer, K. V. Vienola, and J. F. de Boer, “Angiography of the retina and the choroid with phase-resolved oct using interval-optimized backstitched b-scans,” Opt. Express 20, 20516–20534 (2012). 22. K. Nakai, F. Gomi, Y. Ikuno, Y. Yasuno, T. Nouchi, N. Ohguro, and K. Nishida, “Choroidal observations in vogt-koyanagi-harada disease using high-penetration optical coherence tomography.” Graefes Arch. Clin. Exp. Ophthalmol. 250, 1089–1095 (2012). 23. P. Jirarattanasopa, S. Ooto, I. Nakata, A. Tsujikawa, K. Yamashiro, A. Oishi, and N. Yoshimura, “Choroidal thickness, vascular hyperpermeability, and complement factor h in age-related macular degeneration and polypoidal choroidal vasculopathy,” Invest. Ophthalmol. Vis. Sci. 53, 3663–3672 (2012). 24. J.-C. Mwanza, F. E. Sayyad, and D. L. Budenz, “Choroidal thickness in unilateral advanced glaucoma,” Invest. Ophthalmol. Vis. Sci. 53, 6695–6701 (2012). 25. S. Kuroda, Y. Ikuno, Y. Yasuno, K. Nakai, S. Usui, M. Sawa, M. Tsujikawa, F. Gomi, and K. Nishida, “Choroidal thickness in central serous chorioretinopathy.” Retina 33, 302–308 (2013). 26. S. Usui, Y. Ikuno, A. Miki, K. Matsushita, Y. Yasuno, and K. Nishida, “Evaluation of the choroidal thickness using high-penetration optical coherence tomography with long wavelength in highly myopic normal-tension glaucoma.” Am. J. Ophthalmol. 153, 10–6.e1 (2012). 27. J.-C. Mwanza, J. T. Hochberg, M. R. Banitt, W. J. Feuer, and D. L. Budenz, “Lack of association between glaucoma and macular choroidal thickness measured with enhanced depth-imaging optical coherence tomography,” Invest. Ophthalmol. Vis. Sci. 52, 3430–3435 (2011). 28. Z. Q. Yin, Vaegan, T. J. Millar, P. Beaumont, and S. Sarks, “Widespread choroidal insufficiency in primary open-angle glaucoma.” J. Glaucoma 6, 23–32 (1997). 29. M. E. Martı́nez-Pérez, A. D. Hughes, A. V. Stanton, S. A. Thom, and A. A. B. K. H. Parker, “Retinal blood vessel segmentation by means of scale-space analysis and region growing,” Lecture Notes in Comput. Sci. 1679, 90–97 (1999). 30. J. Staal, M. Abramoff, M. Niemeijer, M. Viergever, and B. van Ginneken, “Ridge-based vessel segmentation in color images of the retina,” IEEE Trans. Med. Imag. 23, 501 –509 (2004). 31. W. Cai and A. C. S. Chung, “Multi-resolution vessel segmentation using normalized cuts in retinal images,” Lecture Notes in Comput. Sci. 4191, 928–936 (2006). 32. M. D. Saleh, C. Eswaran, and A. Mueen, “An automated blood vessel segmentation algorithm using histogram equalization and automatic threshold selection,” J. Digit. Imag. 24, 564–572 (20011). 33. E. Ricci and R. Perfetti, “Retinal blood vessel segmentation using line operators and support vector classification.” IEEE Trans. Med. Imag. 26, 1357–1365 (2007). 34. A. S. G. Singh, T. Schmoll, and R. A. Leitgeb, “Segmentation of doppler optical coherence tomography signatures using a support-vector machine.” Biomed. Opt. Express 2, 1328–1339 (2011). 35. L. Zhang, K. Lee, M. Niemeijer, R. F. Mullins, M. Sonka, and M. D. Abràmoff, “Automated segmentation of the choroid from clinical sd-oct,” Invest. Ophthalmol. Vis. Sci. 53, 7510–7519 (2012). 36. V. Kajić, M. Esmaeelpour, C. Glittenberg, M. F. Kraus, J. Honegger, R. Othara, S. Binder, J. G. Fujimoto, #188755 $15.00 USD Received 12 Apr 2013; revised 9 Jun 2013; accepted 10 Jun 2013; published 25 Jun 2013 (C) 2013 OSA 1 July 2013 | Vol. 21, No. 13 | DOI:10.1364/OE.21.015787 | OPTICS EXPRESS 15788 and W. Drexler, “Automated three-dimensional choroidal vessel segmentation of 3d 1060 nm oct retinal data,” Biomed. Opt. Express 4, 134–150 (2013). 37. M. Sohrab, K. Wu, and A. A. Fawzi, “A pilot study of morphometric analysis of choroidal vasculature in vivo using en face optical coherence tomography,” PLoS ONE 7, e48631 (2012). 38. A. N. S. Institute, American National Standard for the Safe Use of Lasers ANSI Z136.1-2007 (American National Standards Institute, New York, 2007). 39. F. C. Crow, “Summed-area tables for texture mapping,” ACM SIGGRAPH Comput. Graphics 18, 207–212 (1984). 40. J. Kittler, J. Illingworth, and J. Föglein, “Threshold selection based on a simple image statistic,” Comput. Vis. Graph. 30, 125–147 (1985). 41. N. Sang, H. Li, W. Peng, and T. Zhang, “Knowledge-based adaptive thresholding segmentation of digital subtraction angiography images,” Imag. Vision Comput. 25, 1263–1270 (2007). 42. N. Otsu, “A threshold selection method from gray-level histograms,” IEEE Trans. Syst., Man, Cybern., Syst. 9, 62–66 (1979). 43. P. Perona and J. Malik, “Scale-space and edge detection using anisotropic diffusion,” IEEE Trans. Pattern Anal. Mach. Intell 629–639 (1990). 44. A. F. Frangi, W. J. Niessen, K. L. Vincken, and M. A. Viergever, “Multiscale vessel enhancement filtering,” Lecture Notes in Computer Science 1496, 130–137 (1998). 45. A. F. Fr


Introduction
The choroid is a highly vascular layer lying between the retinal pigment epithelium (RPE) and the sclera in the posterior eye segment and it accounts for oxygen and nourishment supply to the outer layer of the retina [1,2]. Morphological alteration of the choroidal vasculature plays an important role in the development of some ocular diseases related to circulation abnormalities such as glaucoma and aged-related macular degeneration (AMD) [3][4][5]. Therefore, imaging of the choroidal vasculature has drawn significant interest in ophthalmology.
Indocyanine green angiography (ICGA) is a common method to image the choroidal vasculature [6] and employs near-infrared (NIR) excitation light and NIR-detection for enhanced penetration beyond fluorescein angiography . However, ICGA only provides an en face projection of dye contrast and does not provide good depth resolution. Signal overlapping in different layers will degrade the vessel contrast especially in those areas where the capillaries are dense. Moreover, ICGA requires dye injection, which is uncomfortable for the patient and can cause severe adverse reactions [7,8].
Optical coherence tomography (OCT) has become a powerful imaging modality in ophthalmology by its cross-sectional and three-dimensional (3-D) imaging capability with high resolution [9,10]. Currently clinical ophthalmic OCT employs a probe beam with a 800-nm wavelength band, which is significantly attenuated by the RPE and the choroid because of strong absorption and scattering, and provides only a limited penetration into the choroid. En-hanced depth imaging (EDI-) OCT offers clear images in the choroid by averaging multiple B-scans [11], and it has been utilized to measure the choroidal thickness [11,12]. However, because EDI requires multiple B-scan acquisition at a single retinal position and hence requires relatively long measurement. This nature makes EDI not suitable for 3-D visualization of the choroid.
An alternative method for choroidal imaging is high penetration (HP-) OCT that uses a probe beam at 1-μm wavelength band [13][14][15]. Several extensions of HP-OCT have also been developed to enable detailed and/or quantitative assessment of the choroid, such as choroidal thickness measurement by using birefringence property of sclera [16,17], the choroidal vessel network visualization based on its scattering property [18], and non-invasive choroidal vessel angiography by using Doppler OCT detection [19][20][21].
The structure of the choroid is generally divided in to four layers; Haller's layer, Sattler's layer, Choriocapillaris, and Bruch's membrane from the outer to inner layers [2]. Besides the Bruch's membrane, the other three layers are formed by dense vessels and distinguished from each other by their vessel diameters. Haller's layer consists of blood vessels with large diameters, lying the out most layer of the choroid. Sattler's layer is a layer consisting of blood vessels with medium diameter. Choriocapillaris is immediately posterior to the Bruch's membrane and consists of small capillaries.
Recently, the whole thickness of the choroid has been regarded as an important parameter for diagnosis since it is highly related to several pathological conditions [22][23][24][25][26]. However, it does not fully represent the choroid property and occasionally results in invalidation of diagnosis. For example, it is reported that glaucoma can result in either thinning or thickening of the choroid in different stages [27]. One of the hypotheses explaining this inconsistency is such that different alteration of vasculature occurs at each sub-layer of the choroid [28]. A 3-D visualization and thickness measurement of the choroidal vessels would reveal the detailed variation of the choroid and offer a better understanding of pathology. Hence, segmentation and quantification of the choroidal vessels is critical for investigation and diagnosis of many ocular diseases.
Since vessel segmentation in medical images is an essential step in the diagnosis of various diseases, many algorithms for vessel segmentation of 2-D or 3-D data have been proposed [29,30]. Choroidal vessel segmentation can provide not only an intuitive visualization of the choroidal vasculature but also semi-quantitative morphological information, which is important to increase the understanding of pathogenesis of circulation-related ocular diseases. Generally, vessel segmentation can be achieved through two different strategies; tracking-based methods [31] or window-based methods [32]. Tracking-based methods rely on region tracking with given or detected seeds to extract the vessels, while window-based methods detect the boundary or the ridge of the vessels by using filters. Tracking-based methods are usually time-consuming, while window-based methods often require further refinement such as performance enhancement by using machine-learning algorithms [33,34].
Recently, Li et al. developed an automated choroidal vasculature segmentation method by using seed-detection and region growing approaches and provided 3-D visualization and segmentation of the choroidal vasculature [35]. Kajić et al. presented a method for a fully automated vessel segmentation that employed a multi-scale 3-D edge filter [36]. Although these methods performed well in the literature, the mathematical models and programming of these segmentation algorithms are relatively complicated.
The adaptive threshold method is one of the simplest window-based methods for segmentation. This method classifies the image pixels of a gray-scale image into two classes by using local thresholds determined by the properties of small sub-sections of the image. The backscattering coefficient of blood is much lower than that of other tissues in the choroid, so the intensity of the OCT signal at the vessel region is relatively low. Thus, the adaptive threshold method is a candidate for a simple choroidal vessel segmentation method.
Quantitative analysis of the choroidal vasculature is important for investigation of the pathogenesis and clinical diagnosis of ocular diseases. For example, Sohrab et al. reported a quantitative analysis of the choroidal vasculature in AMD based on the choroidal vessel density that was measured by spectral domain OCT [37]. However, the choroidal vessel features have not been further evaluated.
The measurement of choroidal vessel thickness and discrimination of the internal layers in the choroid are important for the quantitative evaluation of the choroidal vasculature. The vessel diameter distribution can provide valuable information to understand ocular circulation and vasculature. However, this task is challenging because the choroidal vessels are densely packed and are also highly twisted. In this paper, we present an automated framework to segment and quantitatively characterize the choroidal vessels in volumetric OCT obtained by HP-OCT. First, the choroidal vessel is segmented by a window-based vessel segmentation method. This method is a customized multi-scale adaptive thresholding-based segmentation which works on an en face OCT slice at a constant relative depth to the RPE. After the segmentation, the thicknesses of the choroidal vessels are quantified by multi-scale morphological analysis which again works on the en face OCT with a constant depth to the RPE. A generalized schematic is shown in Fig.  1. After these processes, the choroidal vasculature is displayed by both 3-D volume rendering and depth-resolved en face projection. Our proposed approach is applied to normal eyes and visualizes the choroidal vasculature with quantified vessel thickness information. It enables the non-invasive quantitative analysis of choroidal vasculature.

OCT system and measurement protocol
The algorithms described in this paper were designed to work with high-penetration OCT with a probe beam at a 1-μm wavelength band. The details of the system are described elsewhere [19]. In short, the HP-OCT system is a swept-source OCT using a wavelength scanning light source with a center wavelength of 1060 nm, a spectral scanning width of 120 nm, a full-width-halfmaximum bandwidth of 100 nm, and a scanning speed of 100,000 scans/s. Owing to these specifications, the HP-OCT possesses a measurement speed of 100,000 A-lines/s and the -6 dB depth resolution was measured to be 11.0 μm in air (8.0 μm in tissue). The transversal resolution was designed to be 30 μm on the retina.
To obtain an OCT volume, a 6 mm × 6 mm macular region was scanned with 256 B-scan frames. Each frame consists of 2048 A-lines, where the scanning protocol was a horizontal-fast raster scan. These scanning parameters resulted in a voxel size of 2.9 μm (x, horizontal) × 23.4 μm (y, vertical) × 3.7 μm (z, depth) in tissue. The measurement time was 6.6 s/volume with an 80%-duty of a fast saw-tooth horizontal scan.
The probe power on the cornea was around 1.8 mW which is below the ANSI safety standard [38]. The sensitivity was measured to be 99.1 dB in which the recoupling loss at the fiber tip of the sample arm is compensated.

Choroidal vasculature segmentation
With HP-OCT, the choroidal vascular pattern is clearly visualized in the en face slices extracted at the choroid. A customized procedure was developed to segment the 3-D choroidal vasculature. First, the choroid and a part of the sclera are extracted from the OCT volume and the extracted volume is flatten to RPE. So the en face slice of this flatten volume is a slice with a constant relative depth to the RPE and hence is denoted as an RPE-flattened slice. A multi-scale adaptive threshold-based segmentation method was developed and applied to the RPE-flattened slices for choroidal vessel segmentation as well as segmentation artifact removal.

Choroidal volume extraction
To roughly identify the choroidal region and extract the choroidal slices at each depth, Bruch's membrane was segmented automatically. We employed a method similar to that utilized in Ref.
16. This segmentation algorithm segments the pixels with minimum negative gradients along the depth (anterior to the posterior) beneath the RPE in each A-line, where the gradient is defined as the difference in depth of two adjacent pixels in a linear squared-intensity scale. The segmented pixels represent the Bruch's membrane. The region just beneath the Bruch's membrane with a depth size of 370 μm is selected and flatten as the region of interest. Considering a typical maximum choroidal thickness around the macula, the selected volume covers all of the choroid and some of the sclera in most cases. An en face local mean filter with a kernel size of 20 × 20 pixels (58.6 μm (horizontal) × 58.6 μm (depth, in tissue)) is applied to the volume to reduce speckle and to enhance the image quality in en face slices at each depth. Here the local mean filter was utilized rather than widely utilized other sophisticated blurring filters, such as a Gaussian blurring filter. This is because the local mean filtering is significantly faster than other blurring filters especially if one of an optimized implementation method, such as an integral image method [39], is utilized. Figure 2 shows the despeckled slices extracted at different depths in the choroid. The vessel diameters are larger in more posterior slices.
In addition to Bruch's membrane, the anterior edge of the interface between the photoreceptor inner and outer segments (IS/OS) line is segmented by detecting the maximum positive gradient along the depth in the region anterior to Bruch's membrane. By using the locations of the anterior edge of IS/OS and Bruch's membrane, the hyper-reflective complex around the RPE is extracted. These extracted layers will be used later to identify the shadow signal of the retinal vessel as described in Section 3.4.

Adaptive thresholding method
In our method, the blood vessels are segmented by thresholding with a robust automatic threshold determination method [40]. This method is a simple and fast bi-level thresholding of grayscale images and uses customized threshold values for each subregion of the en face image. The threshold ρ(x 0 , y 0 ; w x , w y ) is determined as where x 0 and x i are horizontal pixel indexes and y 0 and y i are vertical pixel indexes.
W (x 0 , y 0 ; w x , w y ) represents the area of the summations with a size of w x μm × w y μm and centered at (x 0 , y 0 ).
where the partial derivatives are numerically obtained with a standard simple derivative kernel consisting of unities and negative unities with a window size of w x /5 and w y /5 for each direction.
The local threshold is calculated at all pixels in each RPE-flattened en face slice. If the pixel intensity is below the threshold, the pixel is classified as a vessel, otherwise it is classified as non-vessel. Finally, we obtain an output binary image in which a pixel value of 1 is assigned to vessel pixels and 0 is for non-vessel pixels. Figure 3(b) shows an example of a binary image obtained by applying the adaptive threshold to the en face OCT slice shown in Fig. 3(a) with a window size of w x = w y = 375 μm (128 × 16 pixels), where the vessel pixels appear in black. Although several blood vessels are segmented in this image, adaptive thresholding cannot properly distinguish the vessel from the background if the local window includes only background. This results in a segmentation error, e.g., as shown in the right middle area of Figure 3(b), where some parts of the non-vessel region are classified as vessel regions. In this particular case, the region without a vasculature pattern is a scleral region. Here we denote these vessel pixels appearing in the non-vessel region as pseudo-vessel pixels.

Busyness filter
In general, at the subregion containing both vessels and the background, the vessel and background regions are quite distinct from each other in the segmented binary image. On the other hand, if the subregion contains only a non-vessel region, the pixel values alternate frequently in space so that the vessel pixels seem more randomly distributed. According to these characteristic appearances of the true and pseudo-vessel pixels, it would be reasonable to use an edge busyness [41], which measures the jittery appearance of edges, to discriminate the non-vessel regions.
In this method, a pixel with a binary value different with at least one of its 8 neighboring pixels are defined as an edge pixel. The busyness of a sub-region, which has the same size with the window of the adaptive thresholding, is then defined as the ratio of the number of edge pixels over the total number of pixels in this sub-region. Figure 3(c) is the busyness distribution of Fig.  3(b). It is obvious that the busyness value is high in the scleral region, i.e., the region without evident vasculature pattern. Thus, the region with segmentation artifacts can be recognized by applying a threshold to the busyness distribution. For example, the red pixels in Fig. 3(d) are the vessel pixels at the region with a busyness greater than 0.2, and are regarded as pseudo-vessel pixels.
In our particular implementation, the vessel-pixels in the region with a busyness greater than 0.2 (for w x = w y = 375 μm), 0.3 (for w x = w y = 188 μm), and 0.4 (for w x = w y = 94 μm) are regarded as pseudo-vessel pixels and their values are replaced with 0 (non-vessel pixel).

Rejection of false vessel pixels
If the window size is too large for the vessel diameter, the contrast variation among the vessels disturbs the segmentation. We empirically found that an adaptive thresholding with a fixed window size cannot properly segment vessels with diameters smaller than 1/5 of the window size. In other words, segmented vessels with diameters smaller than 1/5 of the window size would be artifacts.
A morphological processing-based procedure is then applied to select reliable vessel pixels in the output of the adaptive thresholding with a specific window size of w x = w y = w. At first, small particles which cannot resist two sequential morphological erosion operations with a 3-pixel × 3-pixel structuring element are removed (IMAQ RemoveParticle VI, LabVIEW 2011, National Instruments, TX). The purpose of this operation is to fill the small holes in the large vessels otherwise they can be regarded as two adjacent thin vessels in the following morphological analysis. After this particle removal, a morphological opening operation is performed using a circular structuring element with a diameter of 0.2w. Since the opening operation is erosion followed by dilation, the vessels with a diameter smaller than 0.2w are eliminated, while other vessels are preserved. After the opening operation, vessel particles with Heywood circularity factors lower than 1.5 are removed, where the Heywood circularity factor is a factor defined as the perimeter of the particle divided by the circumference of a circle with the same area as the particle. This process rejects the vessel particles whose shapes are close to round and are not likely vessels. Therefore, with these procedures, misclassified vessel pixels are removed, and true vessel pixels are preserved. By assuming a cylinder with an exactly circular cross-section as a model of blood vessel, the Heywood circularity factor becomes less than 1.5 when the orientation of the blood vessels becomes steeper than approximately 71 degrees against the OCT en face plain. Since choroidal vessels unlikely to have such a steep angle to the OCT en face plain, this circular particle removal process can safely eliminate only the erroneously segmented particles.

Multi-scale choroidal vessel segmentation
The method described in Section 3.2.3 eliminates false vessels, which appear as vessels with a diameter smaller than 1/5 of w. However, it does not necessarily mean that true small vessels are properly segmented. A small window also does not provide proper segmentation for the following reason. The adaptive thresholding method uses the pixel intensity at the edge structure in an image to determine the threshold as suggested by Eq. (1). If the number of true edge pixels that are at the boundary between a blood vessel and background are too few in the window, the adaptive thresholding method cannot provide a proper threshold. Hence, a small window cannot segment blood vessels with large diameters.
Segmentation of blood vessels with a wide range of diameters is necessary because the diameter of a choroidal vessel can range from several micrometers to hundreds of micrometers. This wide diameter-range segmentation is achieved by a multi-scale approach.
In the multi-scale approach, an OCT volume is processed with four window sizes of w = 47 μm, 94 μm, 188 μm, and 375 μm, where w x = w y = w and each window size corresponds to 1/128, 1/64, 1/32, and 1/16 of the transversal imaging field of 6 mm. The four vessel-segmented images are then combined by the following means.
In the combining process, all pixels in the despeckled roughly segmented choroidal volume, namely the output of the process described in Section 3.1, are first classified into two classes of large vessel pixels and middle to small vessel pixels. Because of light penetration in the choroid and the light scattering property of blood, thick choroidal vessels which mainly distribute in Haller's layer appear with lower OCT signal intensity than thin and middle sized vessels which exist close to the RPE. According to this intensity property, we classify pixels as bright pixels, which mainly correspond to the thin and middle sized vessels, or as dark pixels, which mainly correspond to the thick vessels. The classification is done by using an Otsu's threshold [42] that was defined by using all pixels in the roughly segmented choroidal volume at once. The combined vessel-segmented image V c (x i , y i , z i ) is then defined as where (x i , y i , z i ) is the pixel position, I(x i , y i , z i ) is the logarithmic OCT intensity and k * is Otsu's threshold. V (x i , y i , z i ; w) represents the segmented vessel image with a window size of an adaptive threshold of w μm in which the vessel pixels and non-vessel pixels are represented by 1 and 0, respectively. ∪ represents a logical OR operation. Figure 4 shows the vessel segmentation results corresponding to the OCT en face slices shown in Fig. 2.

Retinal vessel shadow rejection
The retinal vessels cast shadows in the choroidal OCT volume. This shadow signal appears in all of the en face slices and is misclassified as a choroidal vessel. Since this misclassification disturbs the quantification of choroidal vasculature (described in Section 4), the retinal vessels should be segmented and removed from the choroidal vasculature. The shadow of the retinal vessels is most clearly visualized in the hyper-reflective complex including the RPE and IS/OS. A projection of the hyper-reflective complex is obtained, as shown in Fig. 5(a), by averaging the intensity in the hyper-reflective complex which has been segmented in Section 3.1. The image is first denoised by an iterative anisotropic diffusion filter [43] defined as where τ is the index of the iteration, where κ is a constant defining the sensitivity to image edges. In our particular case δ = 0.2 and the κ = 0.1 and the diffusion was performed 10 times. After applying the anisotropic diffusion filter, a Frangi filter is applied to further enhance the vessel pattern, where the Frangi filter has provided a likelihood of a region of interest to contain a vessel pattern or other types of ridge patterns by using the eigenvalues of a Hessian matrix [44][45][46]. The nonlinear parameters of the Frangi filter are customized to discriminate the vessels from background by a threshold of 0.5. Thus, the pixel with a Frangi filter output of greater than or equal to 0.5 is labeled as a retinal vessel and value of 1 is assigned. The value of 0 is assigned for the other pixels. The nonlinear parameters of the Frangi filter were adjusted empirically to provide a good segmentation of retinal vessels with a threshold of 0.5. A morphological closing operation with a circular structural element with a diameter of 60 μm is then performed on the binary output to ensure a connective vessel network. Figures 5(b) and 5(c) show the output of the Frangi filter and the binary retinal vessel segmentation result, respectively. The choroidal pixels beneath these segmented retinal vessels are labeled as invalid and excluded from the quantification process described in Section 4. Note that the segmented retinal vessel is slightly more dilated than the real vessel network. This feature further warrants a sufficient removal of the retinal vessel shadow in the choroidal volume.

Methods
In a strict sense, the vessel diameter is determined as the diameter of the blood vessel in a cross-section which is perpendicular to the direction of the vessel. However, the orientation of the choroidal vessel is highly randomly distributed and it is difficult to obtain the cross-section perpendicular to the vessel direction.
In our method, the vessel diameter is determined as the transversal vessel thickness in the segmented en face slices. Similar to what was mentioned in Section 3.2.3, morphological opening was applied to reversed vessel segmented slices in order to eliminate the blood vessels which possess a smaller diameter than the structuring element. Owing to this property, the blood vessels can be classified into two classes by the morphological opening operation; vessels with larger or smaller diameter than the structuring element. According to this property, we can estimate the vessel diameter by applying a series of opening operations with different diameters of circular structuring elements. The diameter of the structuring element in the i-th opening operation is set to be d i = i× d 1 , where d 1 is the structuring element diameter of the first opening operation and the maximum value of d i should be larger than the maximum diameter of the choroidal vessel.
To properly estimate the vessel diameters with different orientation, the structuring element should have a round shape. However, the structuring element is a pixelized binary object and a too small structuring element cannot be approximated as round. To guarantee the round shape of the structuring element with the smallest diameter, the pixel density of the en face slices are increased by 2-D bi-linear interpolation into 2048 pixels × 2048 pixels prior to the opening operation.
After the interpolation, a series of opening operations are applied, where the number of the opening operation was 20, and d 1 = 15 μm. Namely, the diameter of the structuring element varies from 15 μm to 300 μm. Figure 6 shows the diameter estimation results corresponding to the OCT en face slices shown in Fig. 2 where the pixel brightness represents the vessel diameter. Figure 6

A fast approach
In the method described in the previous section (Section 4.1.1), the resolution of the diameter estimation is determined by the diameter increment of the structuring element. Namely, small increments with large number of opening operations provides high diameter-resolution. However, the opening operation is time consuming if the image sizes are large. In the method described in Section 4.1.1, the pixel numbers of both the en face vessel images and the structuring elements are large. Therefore, it often requires a long processing time for vessel diameter quantification. In addition, the processing time for the diameter estimation occupies more than half the entire processing time in our implementation. And hence acceleration of this process is of great interest.
A strategy to accelerate this process is to keep the structuring element small, such that the minimum acceptable size to be regarded as a round shape, and reducing the size of the vessel segmentation image to have a proper relative sizes of the structuring element and the vessel segmentation image. In this fast method, the diameter of the structuring elements is fixed as alternatively, the segmented en face vessel image is resized to N/i × N/i. After applying the opening operation, the resulting image is resized into the original size.
When we set d f = 5 × d 1 and i = 1 · · ·20, the same as the the demonstration in Section 4.1.1, this fast method performed 5-times faster than the standard method in our implementation. Figure 7 shows the difference of two diameter estimation images where one is by the standard method and the other is by the fast method. In most of the regions, the difference is negligible (less than 15 μm). Although the difference becomes relatively large in some thick vessels, it only appears at the edge of the vessels. Since the transversal width of the error is in the same range as the transversal resolution of our OCT, which is around 30 μm, the difference is believed being caused by the limited system resolution, and is not significantly affect to the quantification of the choroidal vasculature.
around 30 μm, which is with a thickness originated from it is in the same range as the The high estimation error at the edge of the large vessels around 30 to 60 μm is also accounted by this mosaic effect. However, this difference would not be significant for the diameter quantification because it is caused by a small difference in the delineated edge of the vasculature between two methods, not by the error in the estimation of the vessel diameter. Since this result has demonstrated sufficient performance of this fast method, we used the fast method in the following sections.

Vascular network thickness
In addition to the choroidal vessel diameter, the thickness of the layer of the choroidal vascular network was also quantified.
Here the vascular thickness is the thickness between the RPE and the posterior envelope of the whole choroidal vasculature, and the envelope is determined by applying an active deformable surface model to the RPE-flattened segmented choroidal vessel volume. We first initialize a deformable plane with 10 × 10 control points as an en face plane as shown in Fig. 8(a). This plane is then dropped from the anterior choroidal region to the posterior direction where the force applied to the each control point is defined as where τ is the iteration count of the position update of the control points and j is the ID of the control point. R defined as where S(x i , y i ) is a smooth surface defined by the control points with bi-cubic interpolation and (u j , v j ) is the transversal position of j-th control points. It should be noted that the positive depth direction represents the posterior direction. P (τ) j is the pushing force defined as the mean vessel diameter in a local region where the local region is a tile centered at the control point and in the same depth with the control point. The size of the local region is 6/9 mm × 6/9 mm, i.e., the same size with a single transversal mesh of the control points. G is a negative constant representing gravity which constantly push the control points into the anterior direction. α and β are predefined positive coefficients to balance these three forces.
The depth positions of the control points are updated for each iteration as where z (τ) j is the depth pixel index of the j-th control point at the τ-th iteration and F θ is a predefined positive threshold of force. The iteration is continued until the depth positions of all control points converge into stable values.
After the convergence, a vascular envelope is obtained by applying 2-D bi-cubic spline interpolation to the control points as shown in Fig. 8(b).

Results
To validate the performance of the algorithms, the macular region of healthy subjects was imaged by HP-OCT with a probe wavelength of 1 μm. The measurement protocol has been described in Section 2.
An example of reconstructed 3-D choroidal vasculature is shown in Fig. 9 (Media 1). The brightness of voxels represents the diameter of the vessel, where thicker vessels appear brighter (a more whitish color). Haller's and Sattler's layers can be distinguished from each other in  Fig. 9(b) such as the thin vessels distributed just beneath the RPE, which represents Sattler's layer, while the thick vessels, i.e. the vessels of Haller's layer, are distributed in the deeper region. Figure 9(d) shows the corresponding vascular network thickness map. Obviously, the entire thickness of the choroidal vascular network layer is highly correlated to the distribution of thick choroidal vessels shown in Fig. 9(a). Figure 10 illustrates the vascular network thickness maps (left) and distributions of the mean choroidal vessel diameter as a function of depth from the RPE in nine different sectors in the maculae (right) of two normal eyes of two subjects. Here the macula was sectorized into nine sectors as a central subfield (C) with 1-mm diameter, nasal inner macula (NI), superior inner macula (SI), temporal inner macula (TI), inferior inner macula (II), nasal outer macula (NO), superior outer macula (NS), temporal outer macula (TI), and inferior outer macula (IO), where the inner and outer rings have 3 mm and 6 mm diameters, respectively. These sectors were originally defined by the Early Treatment Diabetic Study (ETDRS) layout, and have also been used for sectorized choroidal thickness analysis [47].
Obvious increases in vessel diameter along the depth were observed in all of the sections in two eyes. This matches to our anatomical knowledge in the choroid. In case-1, the mean vessel diameter curves show a high gradient at the depth region from 50 to 100 μm. This may be the interface of Sattler's layer and Haller's layer, and the thickness of 50 to 100 μm is consistent with the generally known thickness of Sattler's layer. The thickness of the thick vascular network, which would be Haller's layer, is highly correlated with the vascular network thickness (left images). For instance, it is thick in the SO and IO sectors in Case-1, and thick in the SO and TO sectors and thin in the NO sector in Case-2. It is also noteworthy that the mean choroidal vessel thicknesses in Haller's layer are not so different among the different sectors in the second case.  ICGA provides a two-dimensional visualization of thick choroidal vessels as shown in Fig.  11(a). However, the edges of the thick vessel are unclear because of the blocking effect of anterior thin vessels and dye leakage, and the thin choroidal vessels are hardly visualized because of low contrast. An ICGA-style projection of choroidal vessel diameter volume, denoted as a vessel diameter projection, was created by averaging the diameter value in each A-line. Figure 11(b) shows the vessel diameter projection, in which the pixel brightness is in proportion to the square-root of the depth-averaged diameter value. Since the square-root of the averaged diameter is roughly in proportion to the vessel thickens, it would also be roughly proportional to the amount of dye molecules at each transversal position of the ICGA. And hence the vessel diameter projection would provide a similar image to ICGA. Although similar thick choroidal vessel patterns appeared in the projection and ICGA, the contrast in the projection of the vessel diameter volume clearly possesses better contrast. In the vessel diameter projection, not only are the thick vessels clearly visualized but also thin vessels are seen.
Furthermore, by using the depth-resolved property of OCT, the depth-position can be integrated into the vessel diameter projection as shown in Fig. 11(c). Here the brightness of the pixel represents the vessel thickness the same as Fig. 11(b), while the hue value denotes the diameter-weighted depth position of the vessel defined as where Z d is the diameter-weighted depth position, z i is the pixel index along the depth, M is the number of pixels along the depth, d i (z i ) is the quantified diameter at the depth of z i , and Δz is the pixel separation in micrometers along the depth. The color alteration corresponds to the mean vessel depth weighted by the vessel diameter where red to blue corresponds to from the RPE to the sclera. This en face provides intuitive and simultaneous visualization of the thickness and depth-position of the choroidal vessels.
To qualitatively verify the usefulness of the method, 8 eyes of 4 healthy subjects were examined. The vascular network thickness maps of the 8 eyes are shown in Fig. 12, where each row represents the pair of eyes of a single subject. The vascular network was thinner at the nasal area than the temporal area in all cases. The patterns of the vascular network thickness were found to be symmetric between the right and left eyes except for Subject-3 who had localized thinning at the superio-temporal region in his right eye. Figure 13(a) shows a B-scan corresponding to the dashed line in the vascular network thickness map of Subject-3 in Fig. 12. The region of the localized thinning is indicated by the red circle, and no vessels appear at this region. The same position of the same eye was also examined by polarization sensitive OCT (PS-OCT) [48]. In the phase retardation image taken by PS-OCT ( Fig. 13(b)), increased phase retardation along the depth is observed at the deep region. Since the sclera is highly birefringent [49,50] while the choroid is only moderately birefringent, the chorio-scleral interface is delineated as indicated by the white line. This delineated chorio-scleral interface suggests that the non-vascular region (red arrow in Fig. 12(b)) would not be a part of the sclera but a part of the choroid. This non-vascular choroidal region would partially account for the asymmetric pattern of the vascular network thickness map. Figure 14 shows depth-resolved vascular diameter maps of the 8 eyes of the 4 subjects. It is clear that the vascular network thickness patterns are highly correlated with the thick choroidal vessel distribution.

Discussion
ICGA is commonly used to examine choroid vasculature in the clinic. Although ICGA is safe for most subjects, it is uncomfortable and some adverse reactions limit its application. Although   ICGA is mostly used for visualization of abnormalities in the choroid, it is not suitable for quantitative analysis because of its two-dimensionality and low vessel contrast. Several Doppler OCT-based angiographic methods have visualized 3-D vessel networks in the choroid [19,51]. However, these methods do not visualize the vessels in Haller's layer because of low signal intensity in the vessels and its resulting failure of Doppler measurement. The proposed algorithm in this paper using the scattering OCT obtained by HP-OCT, and hence provides 3-D volumetric visualization of the choroidal vessels as well as all quantitative values of the vessel diameter and the vascular network thickness. Although this method does not provide functional information such as leakage, it would complement other methodologies including FA, ICGA and Doppler OCT.
As discussed in Section 3.3, the choroidal vessel can be properly segmented only if the vessel diameter is around half the window size of the adaptive thresholding. Since the minimum window size in multi-scale adaptive thresholding is 47 μm, the diameter of the most thin vessel that can be properly segmented is around 24 μm. The diameter of the vessel thinner than this minimum diameter would be over-estimated. As a result, the segmentation of the choriocapil-laris is not reliable. Further reductions of the window size cannot improve the accuracy of the segmentation because of the limited transversal resolution of the OCT system. Adaptive optics OCT measurement may supplement the choroidal vasculature analysis.
The precision of vessel diameter quantification is determined by the increment of the structuring element of the opening operation operations. A small increment results in a high precision in vessel diameter estimation, while it results in longer processing time. The trade-off between calculation time and the precision should be balanced to meet the purposes of applications.
The current implementation of the vessel segmentation and quantification was in LabVIEW (LabVIEW 2011 for 64-bit, National Instruments, TX) on a Windows 7 computer with Intel CORE i7 CPU Q720 at 1.6GHz with 8-G RAM. The processing time for vessel segmentation was around 4 s for a single slice of the en face image with 256 × 2048 pixels. Another 20 s was required to quantify the vessel diameter by the standard method using 20 opening operations with initial diameter of 5 pixels and the diameter increment of 5 pixels where the en face image was resized to 2048 × 2048 pixels. This quantification time was reduced to 5 s by using the fast method. In total, the calculation time was 9 s/slice and hence around 20 min/vol where a choroidal volume consists of 100 slices. Although the calculation time is relatively long with the current implementation, the processing time can be accelerated using parallel processing. For example, the adaptive thresholding and morphological opening operation can be accelerated by using a GPU or, alternatively, multi-core CPU processing can parallelize the operations for each en face slice. These parallel processing approaches would reduce the computational time into reasonable ranges for practical clinical use.

Conclusion
We developed an automated method to segment and quantify the choroidal vessels in volumetric scattering tomography obtained by HP-OCT. The three-dimensional choroidal vessel network was segmented based on scattering intensity. The subsequent automated quantification of the vasculature provided both the diameter of each choroidal vessel and the choroidal vascular thickness. Enhanced visualization and characterization of the choroidal vasculature were achieved with normal eyes and demonstrated characteristic patterns of choroidal vasculature.