The retina is regarded as a window to the cardiovascular system. Changes in the retinal microvasculature have been found to be related to several cardiovascular12.–3 and cerebrovascular18.104.22.168.–8 outcomes, among other.910.–11 These evidences make the automatic detection of retinal blood vessels a key step in this area of research. The quantitative description of the detected retinal vasculature can be and has been used to establish the association between retinal vascular properties and clinical and subclinical outcomes, thus providing tools to the clinician for an objective diagnosis.
Extensive work has been done in this field, mainly based on two widely used ocular imaging modalities: color fundus photography (CFP) and fluorescein angiography (FA). Vascular properties such as tortuosity or bifurcation angles computed from these two-dimensional (2-D) fundus images have been associated with several diseases. However, the computation of such properties is incomplete due to the projection to a 2-D plane. A robust method to segment the human retinal vascular network in three-dimensions (3-D) would be a valuable tool and a significant leap forward to fully understanding the pathophysiology of several diseases.
Optical coherence tomography (OCT) is an imaging modality capable of noninvasively imaging the microstructure of tissue in vivo and in situ.12,13 Over time, it became an important tool in the diagnosis of ocular pathologic conditions. It has been extensively used in clinical research and is becoming common in the clinical practice. The principle, based on the backscattering of low-coherence light, is now extensively described in the literature.12,13
On standard spectral-domain OCT (SD-OCT) scans, the retinal blood vessels are not directly visible. Instead, two signatures emerge in the OCT signal. One relates to the fact that hemoglobin absorbs the infrared light. Consequently, backscattering at the structures below perfused vessels is highly attenuated.12,14,15 This effect is well known and has been used to obtain the 2-D vascular network from the 3-D OCT data.1617.–18 While this is a major advantage for 2-D segmentation due to the significant contrast on the retinal pigment epithelium (RPE), the 3-D segmentation of the retinal vasculature requires additional information. The other signature is a diffuse hyper reflectivity on the vessel itself.
In this work, we describe a fully automatic method for the 3-D segmentation of the vascular network of the human retina from standard Cirrus HD-OCT (Carl Zeiss Meditec, Dublin, CA, United States) data and the framework for its reconstruction.
Workflow and Background Work
Following a preliminary study by our research group,22 the 3-D OCT scan is projected to a 2-D ocular fundus reference image. Each pixel in this image translates into an A-scan on the OCT volume. This reference is then used to compute features that are able to discriminate each pixel into vessel or nonvessel groups.23 The A-scans whose pixels were classified into the vessel group, i.e., A-scans that intersect blood vessels, are then processed to determine the depth-wise location of the vessel. Figure 1 shows the global workflow of the algorithm and Fig. 2 shows a graphical depiction of the process.
Throughout this paper, the following coordinate system for the OCT data will be used: is the nasal-temporal direction, is the superior-inferior direction, and is the anterior-posterior direction (depth).
Retinal Layer Segmentation
The very first step consists of determining the depth coordinates of three interfaces: is the inner limiting membrane (ILM), is the junction between the inner and outer photoreceptor segments (IS/OS), and is the interface between the RPE and the choroid. These interfaces are typically the easiest ones to identify.
The retinal layer segmentation process itself will not be discussed here in great detail. On one hand, it does not greatly affect the quality of the results of the final 3-D vascular segmentation and, on the other hand, there is a strong background of published work in this area.24,25 Each B-scan is first filtered with a 2-D Gaussian derivative filter. The rational is that , , and interfaces correspond to the strongest intensity transitions in a normal retina. A threshold is then applied. Finally, 3-D gradient smoothing is applied to the obtained surfaces.
Ocular Fundus Reference Images
As noted, light absorption by hemoglobin is responsible for the decrease in light backscattering beneath perfused vessels. The segmentation process herein takes advantage of this effect by computing a set of 2-D fundus reference images from the 3-D OCT data (Fig. 2). A study was conducted to identify the fundus reference images that provide the best discrimination between the vascular network and the background.26 Four potential reference images were evaluated:26 the mean value fundus (MVF), the expected value fundus (EVF), the error to local median (ELM), and the principal component fundus (PCF).
The PCF image is computed through principal component analysis (PCA) as the principal component of the MVF, EVF, and ELM images. The interfaces at coordinates and are needed for this stage of the process. It was demonstrated that the PCF image provides the greatest extension of the vascular network (equivalent to that achieved with CFP) and the best contrast among the other fundus reference images.26
Feature Computation and Two-Dimensional Classification
To obtain the 2-D vascular network segmentation, we resort to an approach previously published by Rodrigues et al.23 For clarity, we shortly describe in this section the used algorithm and the obtained results.
All four 2-D fundus reference images computed from OCT volumes (MVF, EVF, ELM, and PCF) were used as features in the classification process. Additional features were computed from the PCF image. Specifically, intensity-based features, local-phase features, and features computed with Gaussian-derivative filters, log-Gabor and average filter banks, and bandpass filters were obtained.23 A supervised learning algorithm (support vector machine) was then used to train and classify each A-scan of the OCT volume into vessel or nonvessel A-scans, i.e., to discriminate between A-scans that cross retinal blood vessels and the remaining A-scans. Formally, we define the set of A-scans that cross vessels. The set is defined as the indices of A-scans that cross vessels.
The process proved able to cope with OCTs of both healthy and diseased retinas. It achieved good results for both the macular and optic nerve head regions. For a set of macular OCT scans of healthy subjects and diabetic patients, the algorithm achieved 98% accuracy, 99% specificity, and 83% sensitivity. For both groups, the algorithm favorably compared to the inter-grader agreement.23
Three-Dimensional Segmentation and Reconstruction
The core of the work presented in this paper addresses the issues of estimating vessels diameters and their depth-wise location within the OCT volume, and the 3-D modeling of the vascular network from the human retina, dealing with crossovers and bifurcations. The workflow of this process can be found in Fig. 3.
Vessel Caliber Estimation
The lateral resolution of OCT combined with the optical properties of the vessel walls do not allow for their direct observation. As such, vessel caliber can only be estimated from the shadow cast over the RPE due to light absorption by hemoglobin. Due to the aforementioned reasons (Sec. 2.2), the PCF ocular fundus reference image is the natural choice for the vessel caliber estimation. Thus, it is clear that the estimated caliber is below the real caliber and that this effect is less significant for larger vessels closer to the optic disk, and become more important close to the fovea as vessels become thinner. In addition, the estimated caliber can only be achieved in the plane and we assume them to have a circular cross section.
There is a wealth of literature on methods to compute retinal vessel width from 2-D ocular fundus images. These were tailored for imaging modalities such as CFP and FA.2728.–29 These methods rely on the estimation of the cross section with respect to vessel centerlines and have to deal with the rough definition of vessel borders due to the regular sampling of the digital imaging modalities.30,31 This effect is much more prominent on OCT rendering these methods inappropriate for this modality.
Log-Gabor wavelets are routinely used for vessel enhancement and detection.32 These are created by combining a radial and an angular component (in the frequency space), then determining the scale and the orientation of the filter, respectively. In the time domain, the filter is composed by even (real part) and odd (imaginary part) kernels.
To determine the retinal vessels’ caliber, the PCF image is filtered with a bank of log-Gabor even kernels, each with a unique orientation-scale pair. This bank covers a wide range of scales and orientations to fully describe the whole vascular network. At each vessel pixel , with , the log-Gabor even kernel whose orientation-scale pair better matches the vessel orientation and caliber generates the highest response. The scale of the kernel with the maximal response can then be translated into the caliber of the vessel on that pixel.
The bewildering regions are subsets of neighboring vascular network A-scans where the path of the vessels is unclear from the 2-D segmentation due to bifurcation and/or crossovers (Fig. 4). To determine these regions, several binary morphological operations are used. These regions require further processing (see Sec. 3.4).
First, endpoints in the vicinity of one another on the binary vascular fundus image are bridged and the binary image (the set of vessel points ) is redefined. This updated image is then skeletonized and isolated points are removed. Finally, the bewildering regions are defined as the dilation of branch points and these are removed (erased) from the skeleton image. In consequence, all potential bifurcations and crossovers are adequately left to be re-linked. This new set of points defines the set (Fig. 5).
The set of endpoints on each bewildering region , with , will then be used to look for the most plausible linking configuration of that region (Sec. 3.4).
Vascular Network Depth-Wise Position
As stated, retinal blood vessels are not directly visible in standard OCT data. Typically, vessels on OCT appear as hyper-reflective regions followed posteriorly by the shadowing of structures beneath due to light absorption by hemoglobin in perfused vessels.
In the previous section, all A-scans were classified as vessel or nonvessel A-scans. Only A-scans from the centerline of the vessel are processed to search for the depth-wise location of the vessel, using the following methods.
Principal Component Analysis
PCA is used to enhance similar information between neighbor vessel A-scans in the principal component. For each vessel centerline A-scan (), a circular region in the plane centered at with radius is defined. The set of A-scans in each vessel and nonvessel region is used to compute two new profiles, and , respectively, as follows:
Every A-scan within the defined region is interpolated to account for differences in retinal thickness, from the ILM () to the IS/OS (). One should bear in mind that the whole volume was previously flattened by the IS/OS layer26 and, as such, is now equal for all A-scans.
Each profile in is computed as the principal component of the selected A-scans such that
At the first iteration, the sets are vessel segments delimited by bewildering regions. However, after defining the actual links between segments, the depth-wise position can be recomputed to improve the estimation at the bewildering regions (Fig. 3).
Although the radius was empirically selected following visual inspection of the final 3-D segmentation, it is now kept fixed for the processing of all the examples.
Local Difference Profiles
To estimate the (depth) coordinate of the vessel at each point of the centerline, the difference between the vessel and the nonvessel profiles is computed for each . This operation results in a difference profile with two clear signatures, one due to the presence of the hyper reflectivity and the other due to the shadow in vessel A-scans only, as illustrated in Fig. 6.
Although the vessel walls are not seen in the difference profile, one can estimate their position (and, therefore, the vessel cener) by taking advantage of its diameter estimated in Sec. 3.1. A moving average filter with size allows for the estimation of the center of the vessel by
Bifurcations and Crossovers
Vessel tracking in OCT is quite different from other imaging modalities. Due to the low sampling of OCT systems compared to imaging modalities as CFP and FA, most vessels are just one to two pixels wide and the bewildering regions present many alternatives for the linking process in windows just a few pixels wide (Fig. 4). Furthermore, since we use the shadow to locate the blood vessels, the intensity of different vessels on the fundus reference images do not differ sufficiently to tell them apart.
As a preprocessing step before linkage, we screen all the possible links to find those that are unlikely and to find some links that were not established in the 2-D classification that could help in solving the bewildering region (the so-called growing links). All possible links are then subject to thresholding by length, angle between linked vessel segments, and angle between each vessel segment and the link itself.
For each bewildering region , a cost is assigned to every link between every two points, and , such that and , as
At this point, we shall define all configurations as sets of links subject to the following constraints:
• the set of links contains the link with the least cost ();
• all end-points have to be linked, except points that need vessel segmentation growing to link between them;
• there is only one link by vessel segmentation growing;
• the configuration does not result in intersections or loops in the same vessel.
Generally, the set will enclose all possibilities for crossovers and bifurcations.
The cost of the configuration is then computed as
From the group of many feasible combinations, i.e., ignoring the ones rendering loops within the same vessel, the solution presenting the lowest overall cost is chosen.
In this work, we assume vessels to have a local tubular structure whose centerlines are defined by the 3-D skeleton and the diameter is estimated from the fundus reference image ( plane). In this way, cylinders and cone-like structures are the fundamental components from which the 3-D vascular network is built.
Vessels Path Interpolation
The combination of vessel diameter and low OCT sampling results in vessels is poorly defined in the fundus reference. In addition, at vessel crossovers, the basic assumptions for the determination of the vessel depth location are no longer verified except for the top one. That is, the vessels below the top one do not present the clear signatures (vessel hyper reflectivity and shadow) since they lie on the shadow cast by the vessel above. In this way, interpolation is mandatory and is performed under the assumption that vessels do not present discontinuities or sharp edges, i.e., they are continuous with respect to the first derivative in any of the three dimensions. Under these assumptions, the vascular network is built based on OCT data (control points) making use of cubic spline interpolation.
Vessel reconstruction is achieved by Delaunay triangulation. At each vessel centerline location, the tangent to the path, defined by the centerline points, is computed and the respective normal plane determined. The cross section of the vessel is approximated by a set of points in this cross-sectional plane within the circumference centered at the vessel centerline and a diameter equal to the estimated vessel’s diameter (Fig. 7). Adjacent circumferences are later connected using the Delaunay triangulation in 3-D. The quality of the final reconstruction is directly dependent on the number of triangles, which in turn depends on the circumference sampling and the gap between estimated vessel cross sections.
OCT macular scans of 15 eyes from healthy subjects and eyes from patients diagnosed with type 2 diabetes mellitus (early treatment diabetic retinopathy study levels 10 to 35) were used as the test bed for the proposed methodology. All OCT scans were gathered from our institutional database and were collected by a Cirrus HD-OCT (Carl Zeiss Meditec, Dublin, CA, United States) system. These eyes were also imaged by CFP (Zeiss FF 450 system) and/or FA (Topcon TRC 50DX, Topcon Medical Systems, Inc., Oakland, NJ, United States).
The results obtained for the 3-D vascular segmentation are now presented in three parts: the vessel caliber estimation, the bewildering regions decision, and the reconstruction and position of the vessels.
The automated vessel caliber estimation was compared with the assisted estimates made by three graders (, , and ).
Five vessel segments per OCT eye scan were randomly chosen, each with a minimum of 21 pixels long along the centerline. Graders were instructed to use a software tool. Since the lateral sampling of OCT is relatively low, the graders would not be able to simply mark the vessel boundaries, as is the common practice for CFP.33 Instead, the tool provided for any given caliber, a subpixel continuous marking of the boundaries (for that particular caliber) of the whole vessel segment against the PCF image. Several diameters with increments of 0.2 pixels were tested by the graders to choose the diameter that best fits the data. The result is a set of several caliber measurements (in pixels) for each randomly chosen vessel segment.
For the comparison, two metrics were used, the absolute and the relative differences as defined, respectively, byTable 1).
Comparison between the automatic caliber estimation (S) and each manual grading (G1, G2, and G3), by the mean and standard deviation (SD) values for the absolute difference Da (pixels) and the relative difference Dr between the result G and G^.
Table 1 shows the average metric results by comparing the automatic estimation () to the average grading of all graders and each grader to the average grading of the other two, as a measure of inter-grader agreement. Figure 8 shows the relative difference metric for several ranges of diameters for the same comparisons. As expected, smaller vessels render a higher variability as shown by the histograms. Figure 9 shows both the manual and automatic determination of vessel boudaries over the PCF fundus reference image.
At bifurcation and crossover locations, the continuity of each vessel is hard to detect even for a human grader, naturally depending on the number of possibilities for each region. This was demonstrated to be a very demanding task for the automatic process. For very difficult cases or whenever graders decided, they could access either a fluorescein angiography or CFP, or both, according to the respective availability in our database. In general, the access to these complementary imaging modalities proved to be frequently required for this task. Three metrics were used to assess the accuracy of the automated system in establishing the correct links between vessel segments at bewildering regions. These are: the point accuracy (indicates the percentage of points in correctly grouped), the linking accuracy (the percentage of links correctly connected and nonlinks correctly unconnect, i.e., the percentage of true positives and false negatives) and, the percentage of bewildered regions correctly classified.
Vascular Network Depth-Wise Position
The low visibility of vessel markers (such as the shadow) renders the manual detection of blood vessels in OCT B-scans a hard task. Two graders ( and ) were instructed to mark, at 50 randomly selected vessel A-scans, the position where the shadow of the vessels began. Some restrictions were imposed to the random selection to guarantee an unbiased evaluation: the same vessel could not be selected twice and not more than five A-scans per exam were possible. Furthermore, very small vessels (less than two pixels in radius) were excluded as the graders found it very difficult to evaluate them. Both graders evaluated the same A-scans so that inter-grader variability was computed. A software tool was developed to help in this process. The second grader repeated the process () to establish an intra-grader variability. The results are shown in Fig. 12 and summarized in Table 2.
Comparison between the automatic depth-wise position (S) and each manual grading (G4, G5.1, and G5.2), by the mean and standard deviation (SD) values for the absolute difference Da (mm) between the result G and G^.
The high inter and intra-variability values obtained (0.020 and 0.011 mm, respectively), clearly show how difficult the manual process is.
As is visible from Fig. 12, the automatic estimation of the depth-wise position appears mostly at a higher depth than the position estimated by the graders, since the grader aims at the beginning of the shadow (vessel top wall) and the algorithm aims at the vessel centerline. Hence, this systematic deviation is expected.
The proposed reconstruction seems very feasible. Overall, the reconstructed vessel network is smooth and behaves as expected (Figs. 13 and 14). In fact, the details of the crossovers of vessels, as illustrated in Fig. 11, are also according to physicians’ expectations, where the vessels rapidly deviate to form the crossover. Note that the high axial resolution of OCT is crucial to detect crossovers due to this aspect. Moreover, we illustrate in Fig. 15 the cross section of the vascular reconstruction in several B-scans. These figures illustrate that the method leads to a feasible reconstruction of the position and shape of the vascular network within the OCT volumetric scan.
From Rodrigues et al.,23 the execution time for the 2-D segmentation process (OCT fundus reference computation, features computation, and support vector machines classification), using a MATLAB® (The MathWorks Inc., Natick, MA) implementation, was (). The system hardware used was an Intel® Core i7-3770 CPU (Intel Corporation, Santa Clara, California) at 3.4 GHz.
The additional execution time to achieve the 3-D reconstruction, also using a MATLAB implementation, was () on an Intel® Core i7-4770 CPU at 3.4 GHz. For the 3-D reconstruction, the required time greatly depends on the vessel tree complexity. The high standard deviation value reflects this behavior.
Discussion and Conclusions
The method presents a good overall performance. The location of the vascular tree is (as expected) in the upper third of the retina. The vessel caliber estimation achieves a precision similar to a human grader and the depth position detection is in agreement with the known anatomy. However, the linking of the vessels in the defined bewildering regions works well when few connections are involved, but does not seem to contain enough information to achieve higher accuracy when the number of vessels to connect is higher. Apart form the presented cost functions, many other approaches were tried leading to similar but no better results. Although the bewildering accuracy is about 65%, we note that the linking accuracy and point accuracy are quite high for the problem at hand, bearing in mind that OCT lateral resolution is low (in comparison with other modalities). To the best of our knowledge, this problem is addressed here for the first time for OCT data.
The different steps of the validation show that the algorithm performs well. However, Doppler OCT would be a better ground truth than manual segmentation. Unfortunately, we have no access to such a system at our institutions. The high inter and intra-variability values from the manual gradings indicate how difficult the manual segmentation is.
The first tests on the relevance of the depth component of the vascular network of the retina have already been conducted.33 In this preliminary study, it is shown that the most widely used metric of vessel tortuosity does not have a statistically significant linear relation between 2-D and 3-D metrics. Thus, the use of 3-D tortuosity metrics to correlate with disease status might have a significant impact on the correlation values when compared with those obtained using 2-D metrics.
It is expected that more severe pathologies would render a more difficult task to overcome. In the future, we hope to extend our tests to these cases. Although one might slightly anticipate worse results, please note that our goal of early diagnosis of pathology requires working with eyes within the early stages of disease progression (close to normal).
In the present work, we propose and describe a method to segment and visualize the retinal vascular network in 3-D for OCT data. Increased sampling and accuracy would improve the algorithm’s performance and would allow for an objective validation.
This work was supported by FCT under the research projects (Project Nos. PTDC/SAU-ENB/111139/2009 and PEST-C/SAU/UI3282/2013), and by the COMPETE programs FCOMP-01-0124-FEDER-015712 and FCOMP-01-0124-FEDER-037299. The authors would like to thank António Correia for performing the manual gradings.