Framework for quantitative three-dimensional choroidal vasculature analysis using optical coherence tomography

: Choroidal vasculature plays an important role in the pathogenesis of retinal diseases, such as myopic maculopathy, age-related macular degeneration, diabetic retinopathy, central serous chorioretinopathy, and ocular inflammatory diseases. Current optical coherence tomography (OCT) technology provides three-dimensional visualization of the choroidal angioarchitecture; however, quantitative measures remain challenging. Here, we propose and validate a framework to segment and quantify the choroidal vasculature from a prototype swept-source OCT (PLEX Elite 9000, Carl Zeiss Meditec, USA) using a 3 × 3 mm scan protocol centered on the macula. Enface images referenced from the retinal pigment epithelium were reconstructed from the volumetric data. The boundaries of the choroidal volume were automatically identified by tracking the choroidal vessel feature structure over the depth, and a selective sliding window was applied for segmenting the vessels adaptively from attenuation-corrected enface images. We achieved a segmentation accuracy of 96% ± 1% as compared with manual annotation, and a dice coefficient of 0.83 ± 0.04 for repeatability. Using this framework on both control (0.00 D to − 2.00 D) and highly myopic ( − 8.00 D to − 11.00 D) eyes, we report a decrease in choroidal vessel volume (p < 0.001) in eyes with high myopia.


Introduction
The choroidal vasculature is a rich and complex vascular network that resides between the retinal pigment epithelium (RPE) and the sclera, supplying oxygen and metabolic substance to the photoreceptors and draining away metabolic wastes from the RPE [1,2]. It plays an important role in several ocular diseases such as myopic macular degeneration, age-related macular degeneration (AMD), diabetic retinopathy, central serous chorioretinopathy, and ocular inflammatory diseases. [3,4]. Traditional imaging modalities have been widely used to study choroidal vessels [5], including indocyanine green angiography (ICGA), laser Doppler flowmetry, and laser speckle flowgraphy (LSFG). Unfortunately, these approaches are limited to a two-dimensional (2D) visualization of the projected complex vascular network, and only large choroidal vessels could be dominantly visualized and quantified [6][7][8][9]. Optical coherence tomography (OCT) visualizes the three-dimensional (3D) choroidal vasculature in a dye-free, non-invasive and high-resolution manner [10][11][12]. Given that choroidal vessels reside underneath the RPE, spectral-domain (SD) OCT at the 800nm near Infrared (NIR) window are not ideal for choroidal imaging because it associates with higher scattering loss at RPE. Swept-source (SS) OCT at 1060nm wavelength regime, on the other hand, enables higher signal penetration into the tissue with a lower sensitivity roll-off rate. It yields unprecedented visualization of the choroidal vessels and holds the potential for choroidal vascular network analysis [13].
Segmenting choroidal vessels from an OCT volume is challenging mainly due to three reasons: the anatomic location of the choroid underneath the RPE generally results in low contrast in OCT images; the attenuation differences between the choroidal vessels and stroma results in projection artifacts into the underlying tissue; and the choroidal vessels resemble a complex 3D network with vessel sizes varying significantly in a depth-dependent manner. Many efforts have been devoted to the segmentation of the choroid vessels in OCT. For example, an attenuation correction method has been proposed to remove the projection artifacts cast into the deep tissue, allowing for better delineating the choroidal boundary and choroidal vessels [14,15]. Furthermore, various filters, such as the Frangi filter [16] and Hessian matrix [17], have been tested and applied to enhance the visualization of the vascular structures in OCT images. Nevertheless, these filters are usually not adaptable to features with drastically varying sizes and may not be effective for the segmentation of the entire choroidal vessel network. Alternatively, Duan et al. [18] proposed a multiscale adaptive segmentation algorithm for choroidal vessels from OCT enface images. Each enface image was processed with kernels of four different sizes, and the resultant binarized images were merged adaptively to compute the final vessel mask. More recently, a deep learning-based (RefineNet) choroid vessel segmentation algorithm was reported [19], with an accuracy of 82% compared to the ground truth.
To date, limited studies provide an automatic and reliable volumetric choroidal vascular segmentation with a quantification scheme, which would enable choroidal vascular analysis in large scale clinical studies [20]. Here, we propose a framework to automatically segment and quantify the choroidal vessels using a commercial SS-OCT prototype. A size-selective vessel segmentation algorithm isolated the vessels from the attenuation compensated enface images referenced to RPE. We report the accuracy and repeatability of this segmentation algorithm and propose several metrics derived from the 3D vessel network. We also validated this algorithm in a pilot case-control study comparing the metrics between eyes with high myopia, to those with low to no refractive error (controls).

Study participants and imaging protocol
We recruited 5 male and 10 female subjects (average age: 34.6 ± 14.3 years, range: 20-60 years) from the Singapore National Eye Center using two study protocols (CIRB Ref No.: 2018 and 2019/2245). These protocols were reviewed and approved by the SingHealth Centralized Institutional Review Board and conducted per the Declaration of Helsinki. Written informed consent was obtained from all subjects. The refractive errors in the highly myopic subjects (n = 5) were in the range of −8.00 Diopter (D) to −11.00 D, in contrast to the control subjects (n=10, 0.00D to −2.00D). The volumetric retinal scan data were acquired with a prototype SS-OCT system (PlexElite 9000, Zeiss Meditec, Dublin, CA, US), operating at a sweeping speed of 100 kHz and a central wavelength of 1060 nm with a bandwidth of 100 nm. The axial and lateral

Retinal vessel shadow correction
The shadows of the superficial retinal large vessels could be misclassified as choroidal vessels. To address this issue, we first generated a retinal vessel map by axially averaging a slab located 30-40 µm above the RPE (Fig. 2(a)). The map was then converted into a retinal vessel mask with a Frangi vascular filter [21] and automatic thresholding. Lastly, the missing information in the retinal shadows was restored by locating the best-matched textures from the rest of the image [22,23]. Representative enface choroidal images pre-and post-removal of the shadows of large retinal vessels are shown in Fig. 2(a). It was noted that the shadows of the retinal large vessels were well-compensated (yellow markings in Fig. 2(a)) while maintaining the choroidal vascular structure intact.

Attenuation correction
As light enters the tissue and propagates along with the depth, it attenuates through absorption, scattering, or a combination of both. Therefore, the attenuation properties of the superficial tissue could substantially affect the optical intensity and wavefront reaching the deeper tissue, as well as its appearance in OCT [24]. The heterogeneous attenuation coefficients between stromal and blood vessels result in projection artefacts onto the deeper tissue [25]. We applied a method introduced by Girard et al. [26] to transform the attenuated signal (I ′ z ) into its unattenuated form (I z ) through Eq. (1).
where z is the depth indicator and F is the total depth of the OCT. This transformation was performed along with each A-line scan separately. Since the original images were on a logarithmic scale, normalization and exponentiation were performed ( Fig. 2(b)) before attenuation correction (AC), after which the results were logarithmized and normalized again. In this study, the image was normalized between the minimum and maximum image intensity values. Removing the projection artifact via attenuation correction not only achieved an illuminance-corrected scleral image for better choroid-sclera interface identification (yellow markings in Fig. 2(c)), but also enhanced the contrast of the vessels for more accurate segmentation (red markings in Fig. 2(c)).

Identifying the upper and lower boundaries of choroid on enface images
The choroidal vasculature is sandwiched between the RPE/Bruch's membrane and the sclera. Typically, the upper boundary of the choroid is 10-20 µm below RPE [27], and the lower boundary is referred to as the interface between choroid and sclera. The unique topology of the rich choroidal vasculature enabled us to automatically identify the choroidal boundaries based on the structural similarity measurements [28] -structure similarity index measure (SSIM) between consecutive enface images. As seen in Fig. 2(d), the value of SSIM first decreases as the enface section entered the Bruch's membrane from RPE, and then increased when it entered the choroidal vasculature. In the Haller's layers with larger choroidal vessels, the structural change was slower when traveling along with the depth, reflected by a lowering of SSIM. Approaching the choroidal-scleral interface (CSI), the value of SSIM decreased due to the transition from the vascularized to the non-vascularized region. Within pre-AC images, the value of SSIM stayed relatively consistent due to a persistent shadow appearance of the large choroidal vessels in the sclera (yellow markings in Fig. 2(c)). The AC process enhances this trend, therefore, it facilitates marking the upper and lower boundaries of the choroid vasculature as the intersection points between the two SSIM plots produced by the sequence with and without AC ( Fig. 2(d)). The choroid volume confined between these boundaries was then used for choroid vessel segmentation and its 3D reconstruction.

Choroid vessel segmentation
The choroidal vessels have a wide span of diameter (choriocapillaris < medium-size vessels in Sattler's layer < large-size vessels in Haller's layer), hence we propose a size-selective vessel segmentation strategy. A typical enface image (864 × 864 pixels) is denoted by I M,N , where M and N are the image dimensions in x and y-directions (3000 µm x 3000 µm), respectively. Each enface image was split into k square-shaped sub-images (SI k w ) per the size of the sliding window (w: 125 µm to 1500 µm) and each window had 50% overlapping with its neighbors. The vessels were segmented in each of the sub-images separately ( Fig. 3(a)) using a pixel threshold value (τ k w ) determined by Eq. (2) and (3) [18].
where, e i,j is the gradient of I i,j in x and y-direction. Using τ k w as a threshold, each sub-image (SI k w ) was converted into a binary sub-image (BI k w ). The overlapping region among four sub-images ( Fig. 3(b)) was segmented with τ k w from the sub-image with the highest variance [29]. Sequentially segmenting all the overlapping regions in the enface image resulted in a final binary segmented image (FI w M,N ). For a given enface image, the size of the sliding window was adaptively selected per the vessel size. To estimate the vessel size, the enface images were initially segmented using Otsu's threshold [30]. From the resultant binary image, a Euclidean distance transform was calculated that computes the distance between each of the vessel pixels and the nearest background pixel [31]. Also, the vessel centerlines were extracted by skeletonizing the binary image [32], where the value at each centerline pixel was labeled from the corresponding distance map. In each of the enface images, the maximum (D max ), average (D avg ), and standard deviation (STD) of the vessel size were calculated ( Fig. 4(a)). We assumed that the vessel size monotonically increased along with the depth until it reached its maximum (max(D max )) in the Haller's layer. For frames at depths z ≤ t, we applied a linear regression model to estimate D max p from vessel size STD ( Fig. 4(b) and 4(c)). Here, D max p is the predicted maximum vessel size and D max is the originally measured maximum vessel size. The continuous D max p was divided into 4 groups, and in each group, appropriate window sizes (w) were selected (Eq. (4)). For images at depths z>t, the window size was set to half of the image size (w=1500 µm) to avoid falsely segmenting part of the sclera as vessels (Eq .4). When two window sizes were used, the union of both the sizes were used to achieve the final segmentation CI z M,N . Furthermore, OI z M,N is the segmentation result of the original enface image without sub-image splitting.

3D reconstruction of choroidal vessels
After choroidal vessel segmentation, vessel edges were smoothened by convolving the image with a 3 × 3 pixel square kernel. Thereafter, 3D choroid vessel volume was reconstructed by stacking the segmented enface images. The choroidal vasculature is a network of connected vessels; thus, only the largest interconnected object was kept, and small vessel segments were omitted. Iterative (n=4) 3D convolutions with a voxel of size 3 × 3x3 pixels and a threshold of 0.3 were performed to further smoothen the vessel surface. For 3D visualization and presentation purposes, a volumetric data rendering software (Amira; Thermofisher Scientific, Waltham, MA, USA) was used.

Quantitative metrics
We compared our segmentation with the manually-annotated ground truth using intersection over union (IoU) score, Mathews correlation coefficient (MCC), and accuracy score [33]. The manual annotation was done in MATLAB (academic license version: R2020a by MathWorks, Natick, MA, USA) using an inbuilt interactive polygon tool. Furthermore, the Dice score [17] between two repeated scans was used for repeatability evaluation. Given the anterior eye refractive power induces magnification of the retina [34], image magnification correction was done using the formulation described in [35]. Multiple vascular metrics were extracted, including choroid volume, choroidal vessel volume, vessel diameter, vessel density, and perfusion density. Choroid volume was defined as the outer volume of the choroid that contained the choroid vessels [15]. Analogous to vessel size estimation in an enface image (as described in Section 2.3), vessel diameter was determined by calculating the 3D Euclidean distance transform at each pixel of the skeletonized 3D choroidal vessel network [32]. Vessel and perfusion densities were calculated for both depth-resolved 2D enface images and 3D reconstructed choroidal vessel volume. For any given enface image, enface vessel and perfusion densities were defined as the enface vessel length and enface vessel area per the enface area, respectively. Similarly, volumetric vessel and perfusion densities were defined as the 3D vessel length and 3D vessel volume per the total choroid volume, respectively.

Statistical analysis
For the control and high myopia groups, the data were first tested for normal distribution compliance using the Shapiro-Wilk test. Parametric Student's t-test or non-parametric Mann-Whitney U-test was applied on the normal or non-normal distributed metrics for evaluating their differentiation between these groups. Agreement of vascular metrics was assessed using intraclass correlation coefficients (ICCs), where ICC values less than 0.5, between 0.5 and 0.75, between 0.75 and 0.90, and greater than 0.90 indicate poor, moderate, good, and excellent agreement, respectively [36]. Pearson's correlation analysis was used to measure the statistical linear relationship between metrics among repeated scans [37]. All the data here were reported in a format of mean ± STD. p<0.05 was considered statistically significant. Statistical analysis was performed using an open source Microsoft Excel Data Analysis toolbox (XL Toolbox, version 7.3.4 by Daniel Kraus, Würzburg, Germany). Figure 5(a) compares the choroidal vessel segmentation from an enface image using our proposed methods with previously reported approaches from the literature [18], as well as their agreements with the manually annotated images (n=60, Fig. 5(b)). Otsu's method resulted in over-segmentation in the low contrast regions (red arrows in Fig. 5(c)), especially for small vessels. A local mean-based adaptive threshold method [38] may resolve the small vessels but may over-segment the vessels in the background of the image (red arrows in Fig. 5(d)). Duan et al. [18] used multiple-sized non-overlapping windows to segment the images and jittery segmentation with multiple spurs could be observed (red arrows in Fig. 5(e)). On the contrary, the proposed method adaptively used overlapping windows with different sizes. Figure 5(f-h) demonstrates the effect of window size (w=250µm, 500µm, and 1000µm) on the segmentation. Size-selective windows were needed to segment vessel features with better accuracy (per Eq. (4)): a smaller window size resulted in under-segmentation (red arrows in Fig. 5(f)) and a larger window resulted in over-segmentation (red arrows in Fig. 5(h)). The vessels were well-resolved at different depths throughout the choroid (red arrows in Fig. 5(i)). The proposed method could still suffer from false positives and spurs, which were corrected by morphological operations (red arrows in Fig. 5(j)). Table 1 lists the metrics from images that were manually-annotated or segmented by multiple algorithms. The proposed method achieved the highest accuracy, Mathews correlation coefficient (MCC) and intersection over union (IoU) scores.   [18]. Overlapping window segmentation with square window size, (f) w = 250 µm (g) w = 500 µm (h) w = 1000 µm. Present study segmentation results (i) before and (j) after morphological processing.

Repeatability analysis
We analyzed the repeatability of the proposed segmentation algorithm by evaluating the agreement in depth-resolved enface average vessel size, enface perfusion density, and enface vessel density measured from the repeated scans of three eyes. Figure 6(a-b) shows the 3D segmented choroid vasculature from two repeated scans, visualized from the posterior side. We observed some lateral translations between two scans; however, the overlapped regions showed similar angioarchitectures. Figure 6(c) shows that the depth-resolved mean choroid vessel sizes were highly correlated between the two repetitions. The Bland-Altman plot (Fig. 6(d)) reveals that the average vessel sizes were comparable between two repetitions with >90% of the data points fall within the 1.96 STD range (95% confidence interval). Table 2 summarizes the repeatability test results for three subjects. Excellent repeatability was achieved (ICC >0.94 and R >0.94; p <0.001), and the average volumetric Dice score was found to be 0.84 ± 0.04.

Choroidal vessel network in high myopia eyes
Three-dimensional reconstruction of the choroid vasculature provides a qualitative means to evaluate choroid vessels in the case of high myopia (Visualization 1). As shown in Fig. 7(a) and 7(b), an apparent thinning of the choroid could be observed. This is further corroborated by the quantitative parameters summarized in Fig. 7(c). In contrast to control subjects (Visualization 2), maximum thickness, choroid volume, and choroid vessel volume were significantly lower in high myopia eyes (p<0.05). This is in agreement with the previous literature, where mean sub-foveal choroid thickness decreased in high myopia (93-226 µm) from normal range (203-305 µm) [4,[39][40][41][42]. Also, we found no significant difference in mean and maximum choroidal vessel diameter values between the two groups (p>0.05). Furthermore, in highly myopic eyes, the choroidal volumetric vessel and perfusion densities were found significantly lower (p<0.05) than control eyes.

Discussion
Here, we presented a framework for automatic 3D choroidal vessel network analysis from enface images. This framework included the following steps: removal of retinal vessel shadow from the choroid, correction of attenuation, isolation of the choroid, adaptive segmentation of the choroidal vessels, and reconstruction of the 3D choroidal vasculature. Our framework yielded highly repeatable vascular metrics and we successfully demonstrated its usefulness in a pilot high myopia group.
To date, the segmentation of these "dark vessels" has been based on their negative contrast in OCT intensity images: isolating the dark areas from their surrounding bright areas. With higher penetration and better signal roll-off, swept-source OCT systems could provide better visualization of deeper structures. The resulting improved visualization of the choroidal stroma and sclera will enhance the ability to isolate choroidal vessels. Furthermore, the "darkness" in OCT could also come from shadows underneath the highly attenuating tissue, which can be falsely segmented as choroidal vessels. It is important to include attenuation correction, which has been reported to increase the accuracy of choroid segmentation (91.8 ± 3.7% versus 74.5 ± 8.0%; p< 0.01) [14]. It should also be noted that, although OCT angiography (OCT-A) techniques have been successfully used for retinal vasculature, it was less successful in visualizing vessels with low intrinsic contrast, including lymphatic vessels and choroidal vessels [43]. Particularly in choroidal vessels, OCTA can be used to gain insight into choriocapillaris but not into large vessels [44]. However, choroidal vessels can be imaged with use of exogenous contrast agents or in the absence of RPE due to an increase of vessel SNR in the OCTA images [45].
Selecting an ideal choroidal vessel segmentation scheme in noisy choroidal images is challenging. Using a threshold to distinguish vessels from stroma over the entire database may pose anomalies in accurately segmenting vessels in OCT volume data with uneven contrast and poorly visible vessel regions [46]. Averaging several acquisitions could enhance the choroidal vessel visualization but is consequently affected by inter-scan motion and is more time-consuming [16]. Therefore, a local region-based segmentation approach is more appropriate. This was firstly proposed by Duan et al. [18], where the window size of the local region segmentation was determined by the image depth. Built upon this idea, we further included local threshold segmentation through overlapping windows, vessel size-based adaptive window size usage, and the combination of partial segmentation results from different windows. While the importance of local threshold-based segmentation in the sub-regions of the image has been well-established, the use of overlapping windows ensures continuity in the segmented vessels from an uneven vessel contrast image [29]. However, the downside is the cost of computation time, especially for anterior choroid layers where smaller-sized segmenting windows were adapted. The average computation time for a typical choroidal vessel volume was 18.05 ± 2.07 minutes per choroid volume (100 ± 6 enface images, image size: 864 × 864 pixels), performed on an Intel Core i7-8700K CPU @3.70 GHz system with 64GB RAM. Previously, such computational time was reported to be 20 minutes per volume for a volume consisting of 100 slices (image size: 256 × 2048 pixels) [18]. The processing time can be substantially decreased with the application of GPU-based parallel computation. A deep learning-based segmentation model could be a promising alternative. In the past, a RefineNet model achieved a segmentation accuracy of 82% against manually annotated choroidal vessels in selected B-scans [19]. The downside of the deep learning-based approach is the need for meticulous manual annotation of vessel features, which is, in particular, both challenging and limited for small vessels in the anterior choroid layers.
Validating choroidal vessel segmentation is not a straightforward task. Zhang et al. [17] supported the use of repeatability metrics instead of methods like ICGA, fundus photography, fluorescence angiography, and OCT examination in postmortem eyes to validate choroidal vascular segmentation. Cross modality validation is very challenging, and more importantly, the 2D fundus images could not serve as the ground truth for the 3D OCT scans. Furthermore, choroidal vessel validation against histology is limited by factors such as morphological change due to tissue preparation and a co-registration algorithm design. We validated our segmentation against manual annotation, wherein a segmentation accuracy of 96% was reported. Additionally, we reported an average volumetric Dice coefficient of 0.84 ± 0.04 for the repeatability test, where a previous study reported a Dice coefficient of 0.78 ± 0.08 [17]. It is also worth noting that manual annotation was applied on the attenuation corrected images to avoid potential bias.
We also demonstrated the clinical application of the proposed segmentation framework in myopia. As the role of the choroid in myopia progression and development of pathologic myopia is increasingly being studied [47], accurate segmentation and measurement of the choroid and its layers can be challenging [48,49]. Moreover, while choroidal thinning in myopic eyes has been well documented, this has been limited by quantitative vascular level change analysis [4,[39][40][41][42]50,51]. To the best of our knowledge, only vessel and perfusion densities in 2D enface projection have been reported thus far [52]. In the present pilot study, we have shown that 3D choroidal vessel segmentation could provide multiple volumetric vascular metrics that are significantly affected by high myopia. This demonstrates the potential for studying the topological changes of the choroidal vascular network in myopia and other ocular diseases with choroidal pathology. Our framework that allows for choroidal vasculature quantification could be useful in automatically deriving choroidal measurements in larger cohort studies in the future.
This study has several limitations. First, the implementation is limited to a 3 × 3 mm macular region. This shall be addressed in our future work with the inclusion of a larger field of view data. Second, in myopic eyes, the highly curved surface could distort OCT images, which is not addressed in the current study. However, distortion is minimal in a small FOV [53], and will need to be considered in a larger FOV (12 × 12 mm). Third, we have not tested our algorithm on eyes with pathology in which RPE disruption may be present. Drusen or RPE absence in geographical atrophy can substantially alter the choroidal vessel appearance in OCT images [54][55][56]. Fourth, the current algorithm is not suitable in a real-time clinical setting as it takes approximately 18 minutes to segment a single 3D choroid vessel map. Incorporating GPU application is a promising solution to reduce the computation time; however, this will increase the scan cost. Nonetheless, we have successfully demonstrated this proof of concept with applications to a preliminary clinical study, which could be further developed and improved for clinical application in the future.

Conclusions
In this study, we introduced an automatic 3D choroid vessel segmentation and quantification framework. This framework involved a series of steps: removal of retinal vessel shadows from the choroid, attenuation correction, isolation of the choroid, adaptive segmentation based on vessel size, and 3D reconstruction. High agreement with manual annotation was achieved, along with excellent repeatability. The pilot study using this framework on quantifying choroidal vasculature in high myopia eyes demonstrates its clinical potential.