Age-related macular degeneration (AMD) is the leading cause of irreversible vision loss in older adults in the developed world.1 The neovascular form of AMD is characterized by pathologic new vessels that grow from the choroid, through Bruch’s membrane (BM), and into the avascular outer retina. There, the choroidal neovascularization (CNV) can cause subretinal hemorrhage, fluid exudation, or a number of other deleterious effects that negatively affect vision.2 Fluorescein angiography (FA) and indocyanine green angiography are standard methods used to diagnose and evaluate CNV in the clinic. However, they are two-dimensional and require intravenous dye injection.
Optical coherence tomography (OCT) is a noninvasive imaging technique with micron-scale resolution that is commonly used in ophthalmology to assess retinal morphology. It is based on low-coherence interferometry and produces volumetric images composed of multiple depth-resolved reflectivity profiles. Optical coherence tomography angiography (OCTA) is a functional extension of OCT, which allows for the visualization of vasculature by assessing variation in the OCT signal over time. Regions with high variation represent motion from blood flow while regions with low variation represent static tissue. OCTA has generated significant clinical interest as a label-free method of angiography, and a number of recent studies have investigated its potential in visualizing CNV in AMD.34.5.6.–7 Importantly, OCTA was shown to be able monitor changes in the CNV network over time in response to anti-vascular endothelial growth factor (anti-VEGF) treatment.89.–10 As a result, identification and quantification of the CNV network may have important clinical implications.
Despite rising interests, methods for automated detection and quantification of CNV are lacking. Often, quantification of CNV vessel area is subjective and relies on manual delineation.4,5 To address this issue, we previously developed a saliency-based detection algorithm, which could automatically separate the CNV membrane from noise/artifact on en face angiograms.11 Within the detected CNV membrane, the CNV vessel area was quantified by thresholding to classify pixels as either vessel or noise/artifact and converting the pixel count to CNV vessel area. However, thresholding tended to preserve noise/artifact near vessel edges. More importantly, CNV vessel area gives greater weight to larger abnormal vessels, which may make it a less-sensitive metric to pruning/remodeling of smaller vessel as a result of anti-VEGF treatment.12,13 Thus, we sought to develop a level set method to segment neovascular vessels within the detected CNV membrane and a skeleton algorithm to determine vessel centerlines to quantify the vessel length of the CNV network.
Materials and Methods
Patient Selection and Data Collection
Participants diagnosed with neovascular AMD were recruited from the Casey Eye Institute Retina Service. Diagnosis was based on clinical presentation, a thorough examination, and an FA. Participants were enrolled after informed consent in accordance with an Institutional Review Board/Ethics Committee approved protocol at Oregon Health and Science University and in compliance with the Declaration of Helsinki.
Two volumetric OCTA data sets were collected from the diseased eyes of nine participants using a commercial 70-kHz spectral domain OCT system with a center wavelength of 840 nm (Avanti RTVue XR, Optovue, Fremont, CA). The OCTA scan protocol for a single data set contained two sequential scans, each of which was acquired in less than 3 s and comprised of two repeated cross-sectional B-scans at 304 locations covering a area with a depth of 1.6 mm. Each B-scan was made up of 304 A-scans; each A-scan was 512 pixels in depth. The split-spectrum amplitude-decorrelation angiography (SSADA) algorithm was used to detect flow.14,15 Briefly, SSADA contrasts blood flow from static tissues by assessing the reflectance amplitude decorrelation on repeated B-scans at the same location. Split-spectrum processing was used to improve the signal-to-noise ratio of flow detection. For the first scan in the OCTA scan protocol, the fast scanning direction was in the transverse -direction. For the second scan, the fast scanning direction was in the transverse -direction. The two orthogonal scans were then registered and merged to correct for motion from microsaccades and form a single data set.16
The flowchart of the proposed algorithm is shown in Fig. 1. Following semiautomated segmentation to produce en face angiograms,17 the detection and quantification of the CNV network was achieved in four steps. First, the previously described saliency-based method was used to detect the CNV membrane. Second, thresholding was used to identify the initial CNV vascular region. Next, the CNV vasculature was refined using a level set method. Finally, a skeleton algorithm was used to locate the vessel centerlines. Sections 2.3–2.5 will describe the steps in more detail. The algorithm was implemented with custom software written in MATLAB® 2011a (MathWorks, Natick, MA).
Determining the Initial CNV Vascular Region
To visualize inner and outer retinal circulations, en face angiograms were constructed by maximum flow projection within slabs defined by segmented boundaries using semiautomated software.17 The inner retinal slab includes tissue layers between the inner limiting membrane and outer boundary of the outer plexiform layer (OPL), and the outer retinal slab includes layers between the outer boundary of OPL and BM.
In healthy eyes, the outer retina is typically avascular. However, OCTA suffers from shadowgraphic flow projection artifacts, whereby flow from superficial circulations can be replicated on deeper structures. Application of our previously developed saliency-based algorithm allowed for detection of the CNV membrane in the en face outer retinal angiogram.11 Thresholding was then used to try and separate vessels from noise/artifact to generate the initial CNV vascular region. The threshold was set at two standard deviations above the mean decorrelation value in the foveal avascular zone.
Vascular Boundary Delineation
Thresholding tended to preserve noise/artifact near vessel edges. The flow signal within the initial CNV vascular region was thus inhomogeneous. We used a level set method with a local intensity clustering criterion to better separate the flow signal belonging to vessels from remaining noise/artifact.18 Image inhomogeneity could be described as18
The above-defined local clustering criterion function classified the flow signal only within the circular neighborhood . The integral of the local clustering criterion function for all over the domain constituted the overall function . The minimization of meant the joint minimization of for all over the domain and was done using a level set approach. In the level set method, was the level set function. The vessel region was and the background region . Their corresponding membership functions were defined as , where and . was a smoothed Heaviside function.19 The level set formulation was then as follows:18
Two regularization terms were combined with to generate the variational level set formulation as follows:
The second term on the right side of the equation computed the length of the zero level set function and was used to smooth the zero level set contour.19 The third term was the distance regularization term, and it was used to maintain the signed distance property during the level set evolution.20,21
Pixel classification into or was achieved by minimizing the function in an iterative process. A standard gradient descent approach based on the Gâteaux derivative was applied to minimize the function with respect to the variable for fixed and . The gradient flow formulation was18
In our implementation, the truncated Gaussian function had a scale parameter , and the mask size of the convolution kernel was . The radius was then 8.5. Other parameters were as follows: in the smoothed Heaviside function , the coefficient which modifies the first regularization term in Eq. (5) , which is a scale parameter that modifies the region size for the second regularization term in Eq. (5), and for the gradient flow step size.
Vessel Skeleton and Quantification
The level set method separated vessels from background noise/artifacts to produce the CNV vascular region. This then allowed for the application of a skeleton algorithm on the binary image of the CNV vasculature. We used the Zhang–Suen22 parallel thinning algorithm and implemented it over a window. Briefly, the algorithm iteratively removes pixels until only a skeleton remains. In the first subiteration, south–east boundary and north–west corner pixels are considered. In the second subiteration, north–west boundary and south–east corner pixels are considered.
Quantification of CNV vessel area and vessel length was then performed on the binary image of the CNV vascular region and CNV skeleton, respectively. CNV vessel area was the sum of nonzero pixels converted to . Vessel length was the sum of the pixels making up the skeleton converted to mm.
Verification of Results
The automated algorithm was applied to both OCTA scans from an eye, and within-visit repeatability was assessed using coefficient of variation. The CNV skeleton from the algorithm was compared to manual delineation of the CNV centerlines done by an experienced grader. The expert grader delineated the vessel centerlines of the CNV from the en face angiogram of the outer retinal slab. The grader had access to the inner retinal angiogram for reference. The Jaccard similarity metric was used to compare results from the automated algorithm and manual delineation:
The en face outer retinal angiogram from a participant with neovascular AMD was used as an example to show the workflow of the algorithm. Figure 2(a) shows the original outer retinal angiogram with CNV and projection artifacts. The CNV membrane [Fig. 2(b)] could be detected and extracted using a saliency-based approach.11 We then used thresholding within the CNV membrane to generate the initial CNV vascular region shown [Fig. 2(c)]. The CNV vascular region contained noise/artifacts which would result in overestimation when quantifying CNV vessel area. The binary image of the CNV and its boundaries [red contours in Fig. 2(c)] were used to initialize the level set method. The level set method iteratively performs the following: (1) estimate the flow signal in the vessel () and background () regions based on Eq. (7) (where in the first iteration is the image used to initialize the level set method); (2) reclassify pixels as belonging to vessel or background based on Eq. (6); and (3) estimate the bias field () based on Eq. (8). After the level set function converged ( and change becomes minimal), the final vessel boundaries were smooth and not surrounded by obvious noise/artifact [Fig. 2(d)]. The CNV vasculature shown in Fig. 2(e) was binarized and skeletonized. Figure 2(f) shows the vessel skeletons used to analyze the vessel length. The level set method and skeleton algorithm required 13 s to execute on an Intel Xeon CPU (E3-1226, 3.3 GHz); 93% of that time was spent on the level set evolution.
Data from nine participants with neovascular AMD were analyzed. Two volumetric data sets from each eye were evaluated to assess within-visit repeatability. Figures 3 and 4 show the results from one volumetric data set from the other eight participants. The cases included both type I and type II CNV with a wide range of sizes. Compared to determining CNV vessel area by thresholding, the level set method reduced the CNV vessel area measurement by ( deviation). The level set method improved within-visit repeatability, from a coefficient of variation of 7.3% with thresholding to 5.9%.
The skeleton algorithm was then applied to the initial CNV vascular region (thresholding vessel skeleton) and the output after level set evolution (level set vessel skeleton). The vessel skeleton results were compared to manual delineation of vessel centerlines. The level set vessel skeletons showed good agreement with manual delineation, with a Jaccard similarity metric of . This was an improvement (paired -test ) over thresholding vessel skeletons compared to manual delineation, which had a Jaccard of . A comparison of false-positive and false-negative errors (Table 1) showed improvement in both with the level set approach (). A greater change was seen in the false-positive error metric, indicating that the level set method removed noise/artifact pixels to prevent an overly complex vessel skeleton. The vessel length within-visit repeatability as assessed with coefficient of variation was 1.8% for the level set vessel skeleton and 8.1% for manual delineation.
Agreement between automated algorithms and manual delineation of neovascularisation.
|Jaccard similarity metric|
Mean±standard deviation of the Jaccard similarity metric and error rates were computed from the nine participants.
OCTA is a relatively new approach to assessing CNV, and quantitative analysis of the CNV network may have important clinical implications. In this work, we developed a level set method to segment neovascular vessels and skeleton algorithm to determine vessel centerlines. The algorithm first identifies the CNV membrane using saliency-based detection. Next, thresholding is used to obtain the initial CNV vascular region. A level set method then refines the CNV vasculature, using local information to remove remaining noise/artifacts. Finally, a skeleton algorithm is applied to delineate the vessel centerlines based on the output from the level set method. The results of the algorithm as applied to nine eyes with neovascular AMD were shown. A comparison of the vessel skeleton from the algorithm to manual delineation showed good agreement with an average Jaccard similarity metric of 0.69.
Compared to thresholding when calculating CNV vessel area, the level set method reduced the value of the measurement by 24.6%. Inspection of the images showed that this reduction generally was removal of noise/artifact pixels near the edge of vessels. Unfortunately, it also removed pixels that may have been vessels with weak flow. For example (Fig. 5), in the participant shown in Fig. 2, there was removal of noise pixels near the edge of the CNV [green arrows in Figs. 5(b) and 5(c)], but also removal of vascular pixels within loops of the CNV network [yellow arrows in Figs. 5(b) and 5(c)]. A more quantitative evaluation of the impact of the level set method on CNV vessel area may prove insightful for future adjustments or modifications to the level set algorithm. Nevertheless, the application of the level set method to refine the CNV vasculature region simplified and improved the output from the skeleton algorithm [compare Figs. 5(d)–5(f)].
The CNV skeleton derived from skeletonization of the output from the level set method was generally comparable to manual delineation. Discrepancies occurred when the CNV detection missed parts of the CNV that had weak flow signal. For example, in participant 6, the grader drew vessels near the bottom of the CNV network [Figs. 6(a) and 6(b)] that were not included in the initial CNV vascular region [Fig. 4(b2)]. In addition, how the grader interpreted the vessel centerlines within a mass of flow signal did not always align with the algorithm output [Figs. 6(c) and 6(d)]. These issues suggest avenues for future improvement.
While the algorithm was capable of quantifying CNV vessel area and vessel length in the participants of this small study, additional studies with a larger number of diseased eyes are needed. Such a study would provide a better evaluation of the repeatability of the algorithm. Furthermore, it would be particularly insightful to test the algorithm’s ability to detect changes in the CNV network as part of the natural progression of the disease or, alternatively, in response to treatment.
We developed an automated algorithm based on level set segmentation and skeletonization to delineate centerlines of CNV as visualized on en face OCT angiograms. Level set segmentation to refine CNV vessel edges preceding skeletonization was essential to improving the vessel skeleton output. CNV vessel skeleton results from the algorithm applied to OCTA scans from participants with neovascular AMD showed good within-visit repeatability and agreement with manual delineation.
This work was supported by NIH Grant Nos. R01 EY023285, R01 EY024544, DP3 DK104397, P30 EY010572; the National Natural Science Foundation of China (Grant No. 61471226); the Natural Science Foundation for Distinguished Young Scholars of Shandong Province (Grant No. JQ201516); and an unrestricted grant from Research to Prevent Blindness. Oregon Health and Science University, Yali Jia, and David Huang have a significant financial interest in Optovue. David Huang also has a financial interest in Carl Zeiss Meditec. These potential conflicts of interest have been reviewed and managed by Oregon Health & Science University. Other authors do not have financial interest in the subject of this article.
Simon S. Gao received his PhD in bioengineering from Rice University in 2014. He is a postdoctoral fellow at the Center for Ophthalmic Optics and Lasers Laboratory, Oregon Health and Science University (OHSU). His research interests include OCT, OCT angiography, and image processing.
Li Liu received her Master's degree from Shandong University. Currently, she is a PhD candidate in Southeast University. Her research interests include medical image processing and its clinical application.
Dengwang Li received his PhD from Shandong University through a joint PhD program with the University of Sydney sponsored by CSC. Currently, he is an associate professor in Shandong Normal University and is the vice director in the Shandong Province Key Laboratory of Medical Physics and Image Processing Technology. His research interests include medical image processing and its clinical applications.
Yali Jia received her PhD in biomedical engineering from OHSU in 2010. She is an assistant professor of ophthalmology at OHSU. Her research interests include OCT and OCT angiography as well as their applications in retinal diseases.