1.

Introduction

The automatic detection and quantitative description of retinal blood vessels is an important area of research. It can assist the clinician toward the objective diagnosis of vascular pathologies such as retinopathy of prematurity^1,2 or hypertensive retinopathy.³

The retinal vascular network is frequently used as the anchor to coregister data between different imaging modalities.⁴ Moreover, the coregistration of exams taken at different times allows a more accurate study of the development and progression of the pathological process and a deeper insight of the changes occurring in the human retina.⁵

Special focus has been given to the imaging modalities of color fundus photography (CFP) and fluorescein angiography (FA). Tracking-based methods,^6,7 matched filters,^8–12 Hessian matrix and gradient vector fields,^13–15 supervised learning,^11,16–18 and other strategies¹⁹ are several of the many approaches proposed for segmentation of the vascular network from these imaging modalities.

The working principle of optical coherence tomography (OCT), which is based on the backscattering of a low-coherence light, has been extensively described in the literature.^20–22 OCT has made it possible to acquire three-dimensional (3-D) data of the microstructure of biological tissue in vivo and in situ. Its main application is imaging the human retina, for which it has become an important tool. With the introduction of high-definition spectral-domain OCT, it is possible to acquire high-resolution cross-sectional scans while maintaining the acquisition time.

Vessel segmentation from OCT data is considerably different than that from CFP or FA. At the wavelengths used by regular OCT systems (800 to 900 nm), light is absorbed by the hemoglobin, leading to a decrease of the backscattered light from the structures beneath the perfused vessels.²³ As such, segmentation techniques rely on this well-known effect.

A two-dimensional (2-D) image can be obtained by projecting the OCT volume depth-wise. However, traditional methods translate into a vascular network with suboptimal detail and contrast. The lateral resolution and spatial sampling interval of OCT make the process of vascular network segmentation difficult, particularly at the macula.²⁴ Because the vessels are considerably thinner in this region, they present lower levels of contrast.

A robust method for vascular segmentation on 2-D ocular fundus reference images obtained from OCT data would be valuable and a significant starting point toward several algorithms, such as image coregistration and 3-D vascular segmentation algorithms.²⁵

Only a few algorithms for segmenting the retinal vascular network using OCT have been described in the literature. The first method was proposed by Niemeijer et al.²⁴ for optic nerve head (ONH) OCT scans. A 2-D projection (by depth-wise sum) of data from the retinal pigment epithelium (RPE) was used, followed by a supervised pixel classification. Xu et al.²⁶ presented a method that does not require OCT layer/interface segmentation. A boosting learning algorithm was used to segment the vascular network from ONH OCT volumes. Pilch et al.²⁷ took a different approach and segmented the vessels directly on high-resolution cross-sectional B-scans close to the ONH. All of these techniques were validated on healthy retinas.

This report describes a fully automatic method for segmenting the vascular network of the human retina from standard OCT data. We rely on work previously developed by our research group²⁸ to generate 2-D fundus reference images from OCT data volumes. From these images, a set of features are computed to feed a supervised classification algorithm, support vector machine (SVM), to classify pixels into vessel or nonvessel.

2.

2-D Fundus Images

As noted, light absorption by hemoglobin is responsible for the decrease in light scattering 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 at the preprocessing step. A study was conducted to identify the images that provide the best discrimination between the vascular network and the remaining tissue (background). The images evaluated in this study were the mean value fundus (MVF), the expected value fundus (EVF), the error to local median (ELM), and the principal component fundus (PCF).²⁸

Throughout this paper, we use the following coordinate system for the OCT data: $x$ is the nasal-temporal direction, $y$ is the superior-inferior direction, and $z$ is the anterior-posterior direction (depth).

Let $V$ be a low-pass filtered OCT volume, flattened at the junction of inner and outer photoreceptor segments (IS/OS). The MVF image is computed as the average of the A-scan values within the lower layers of the retina, i.e.,

The EVF image of order $\rho$ is computed by

Both MVF and EVF are corrected for nonuniform intensities.²⁸

The ELM image of order $\tau$ is given by

The PCF image is computed as the principal component (by principal component analysis) of the MVF, EVF, and ELM images. In Ref. 28, it was demonstrated that the PCF image provides a greater extension of the vascular network and better contrast than the other fundus reference images (MVF, EVF, and ELM). In addition, when computed from standard OCT data, it presents a vascular network extension similar to that achieved with CFP.²⁸

Figure 1 shows the four fundus images from the same OCT scan covering the central 20 deg field of view of a healthy retina.

Fig. 1

Two-dimensional fundus reference images computed from the same optical coherence tomography ocular scan: (a) mean value fundus, (b) expected value fundus, (c) error to local median, and (d) principal component fundus.

3.

Features

All four 2-D fundus references images computed from OCT volumes (MVF, EVF, ELM, and PCF) were used as features in the classification process. This section describes additional features that were computed from the PCF image. These features were selected through a forward-selection approach from a larger pool of features (see Sec. 5).

Sliding windows and filters discussed in this section consider the differences in spatial sampling between the $x$ and $y$ directions. As such, the results are independent of the acquisition protocol used.

The parameters used for computing the fundus images and SVM features presented herein are discussed in Sec. 7. Throughout the following sections, $n$ and $o$ represent the scale and orientation indexes, respectively.

3.1.

Intensity-Based Features

Intensity-based features (previously used in a similar context¹⁶) are used to describe the local intensity variations. These features are computed from the PCF image $I$ using a sliding window $W(x,y)$ of size $w[\begin{array}{cc}{N}_x& N_y\end{array}]$ and centered at $(x,y)$. The range, average, standard deviation, and entropy features (Fig. 2) are computed as follows:

Eq. (4)

$$I^{\text{range}}(x,y)=I(x,y)-[\underset{k,l\in W(x,y)}{\max }I(k,l)-\underset{k,l\in W(x,y)}{\min }I(k,l)],$$

Eq. (5)

$$I^{\mathrm{avg}}(x,y)=I(x,y)-\frac{1}{N_W}\sum _{k,l\in W(x,y)}I(k,l),$$

Eq. (6)

$$I^{\mathrm{std}}(x,y)=\sqrt{\frac{1}{N_W}\sum _{k,l\in W(x,y)}{[I^{\mathrm{avg}}(k,l)]}^2},$$

Eq. (7)

$$I^{\text{entropy}}(x,y)=-\sum _{k,l\in W(x,y)}I(k,l){\log }_{2}I(k,l),$$

where $N_W$ is the total number of elements in $W$ ($N_W=w^2{N}_x{N}_y$).

Fig. 2

Intensity-based features computed from the principal component fundus image: (a) range, (b) average, (c) standard deviation, and (d) entropy.

3.2.

Gaussian-Derivative Filters

The use of the Hessian matrix to perform scale-space analysis is a well-established technique that is commonly applied in CFP and FA image analysis.^13,14 It is applied here, with slight modifications, to the PCF image.

The Hessian matrix $H$ at scale $n$ is a square matrix defined by

Eq. (8)

$$H_n(x,y)=\left[\begin{array}{cc}{L}_n^{xx}(x,y)& L_n^{xy}(x,y)\\ L_n^{yx}(x,y)& L_n^{yy}(x,y)\end{array}\right],$$

where $L_n^{xx}$, $L_n^{xy}$, $L_n^{yx}$, and $L_n^{yy}$ are the second-order partial derivatives at scale $n$. These derivatives result from the convolution of second-order partial derivatives of a Gaussian filter ($G_n$) with the fundus reference $I$, leading to

Eq. (9)

$$\begin{gathered} L_n^{xx}(x,y)=I(x,y)*{\sigma }_{n_x}{\sigma }_{n_y}\frac{{\partial }^2}{\partial x^2}{G}_n(x,y) \\ {L}_n^{yy}(x,y)=I(x,y)*{\sigma }_{n_x}{\sigma }_{n_y}\frac{{\partial }^2}{\partial y^2}{G}_n(x,y) \\ {L}_n^{xy}(x,y)=L_n^{yx}(x,y)=I(x,y)*{\sigma }_{n_x}{\sigma }_{n_y}\frac{{\partial }^2}{\partial x\partial y}{G}_n(x,y), \end{gathered}$$

where * is the convolution operator and $G_n$ is defined as

Eq. (10)

$$G_n(x,y)=\frac{1}{2\pi {\sigma }_{n_x}{\sigma }_{n_y}}\exp [-\frac{1}{2}(\frac{x^2}{{\sigma }_{n_x}^2}+\frac{y^2}{{\sigma }_{n_y}^2}\left)\right].$$

To accommodate for differences in sampling, the standard deviations are defined as

Eq. (11)

$${\sigma }_{n_j}={\sigma }_n\alpha p_j,$$

where $p_j$ is the spatial sampling along the $j=\{x,y\}$ direction, $\alpha$ is a normalization parameter, and ${\sigma }_n$ is the sampling-invariant standard deviation of the filter at scale $n$.

3.2.1.

Hessian eigenvectors

Eigenvectors from the Hessian matrix $H_n(x,y)$ provide information about the image curvature at $(x,y)$. For ocular fundus images, at vessel pixels, the eigenvector $v=(\lambda ,\theta )$ (in polar coordinates) associated with the largest/smallest eigenvalue (${\lambda }^{+}$ and ${\lambda }^{-}$, respectively) is normal/parallel to the vessel.¹⁴ In addition, ${\lambda }^{+}$ at vessel pixels is considerably larger than the ones in the background, and it presents a consistent direction ${\theta }^{+}$ across different scales. As such, relevant local information can be extracted by resorting to eigenvectors from the Hessian matrix. In this work, we use ${\lambda }^{+}$ (eigenvalue) and ${\theta }^{+}$ (eigenvector orientation—mapped to the interval $]-\pi /2,\pi /2]$) to compute two features.

3.2.2.

Laplacian of Gaussian

The Laplacian of Gaussian filter is commonly used as an edge detector. This feature represents the Laplacian of a low-pass (Gaussian) filtered version of image $I$; it can be directly obtained as the trace of the Hessian matrix for multiple scales

Eq. (12)

$$\mathrm{\Delta }{I}_n=\mathrm{\Delta }(I*G_n)=L_n^{xx}+L_n^{yy}\mathrm{.}$$

3.2.3.

Gradient norm

The norm of the image gradient is computed from a Gaussian low-pass filtered version of image $I$, for multiple scales, as

Eq. (13)

$${|\nabla I|}_n=\sqrt{{(L_n^x)}^2+{(L_n^y)}^2},$$

with

Eq. (14)

$$\begin{gathered} L_n^x(x,y)=I(x,y)*{\sigma }_{n_x}\frac{\partial }{\partial x}{G}_n(x,y) \\ {L}_n^y(x,y)=I(x,y)*{\sigma }_{n_y}\frac{\partial }{\partial y}{G}_n(x,y). \end{gathered}$$

3.2.4.

Features across acales

To reduce the number of features and create new features that encode multiscale information, the results across scales are combined¹⁴

Eq. (15)

$${\lambda }^{\max }=\underset{n}{\max }\left(\frac{{\lambda }_n^{+}}{{\sigma }_n}\right)$$

Eq. (16)

$${\theta }^{\mathrm{var}}=\underset{n}{\mathrm{var}}{\theta }_n^{+},$$

Eq. (17)

$$\mathrm{\Delta }{I}^{\max }=\underset{n}{\max }\left(\frac{\mathrm{\Delta }{I}_n}{{\sigma }_n}\right),$$

Eq. (18)

$${|\nabla I|}^{\max }=\underset{n}{\max }\left(\frac{{|\nabla I|}_n}{{\sigma }_n}\right).$$

Results can be seen in Fig. 3.

Fig. 3

Gaussian-derivative features computed from the principal component fundus image: (a) Hessian eigenvalue, (b) Hessian eigenvector orientation, (c) Laplacian of Gaussian, and (d) gradient norm.

3.3.

Local-Phase Features

Good results in edge and corner detection using phase congruency were achieved by Kovesi.^29,30 The method proved to be useful for vascular segmentation in ocular fundus images.³¹

Computing local-phase features requires the convolution of the image with a bank of log-Gabor kernels, each having a unique orientation-scale pair. These filters are created by combining a radial and an angular component, limiting the frequency bands and the orientation of the filter, respectively (Fig. 4).

Fig. 4

Log-Gabor filter (frequency domain): (a) radial and (b) angular components.

The radial component is computed with the log-Gabor transfer function in the frequency space $F(\varrho ,\theta )$ (in polar coordinates) as

Eq. (19)

$$F_n^{\mathrm{lG}}(\varrho ,\theta )=\exp \left[\frac{-{\ln }^2(\varrho /f_n^0)}{2{\ln }^2\varsigma }\right],$$

where $\varsigma$ is the ratio between the standard deviation of the Gaussian describing the log-Gabor transfer function (in the frequency domain) and the central frequency $f_n^0$, given by

Eq. (20)

$$f_n^0=f^{\max }{2.1}^{-(n-1)},$$

where $f^{\max }$ is the largest central frequency.³² The log-Gabor transfer function is then multiplied (Hadamard product) with a low-pass filter (Butterworth) to ensure uniform coverage in all orientations

Eq. (21)

$$F^{\mathrm{lp}}(\varrho ,\theta )={[1+{(\varrho /f^c)}^{2\upsilon }]}^{-1},$$

where $f^c\in [0,0.5]$ is the cutoff frequency and $\upsilon$ is the order of the filter.

In turn, the angular component is given by

Eq. (22)

$$F_o^{\mathrm{as}}(\varrho ,\theta )=\{\begin{array}{ll}(1+\cos d)/2& \text{if}d\le \pi \\ 0& \text{if}d>\pi \end{array},$$

with

Eq. (23)

$$d=|\theta -{\theta }_o|N_o/2,$$

where ${\theta }_o$ is the orientation and $N_o$ is the total number of orientations.

These components are multiplied (Hadamard product) to obtain the frequency domain log-Gabor filter. In the time domain, the filter is composed by the even (real part) and the odd (imaginary part) kernels (Fig. 5).

Fig. 5

Log-Gabor filter (time domain) for a specific orientation: (a) even and (b) odd kernels.

The local-phase features—phase congruency (PC), , feature type (FT), phase symmetry (PSym), and symmetry energy (SymE)—are described below Examples of these features are shown in Fig. 6.

Fig. 6

Local-phase features computed from the principal component fundus image: (a) phase congruency, (b) feature type, (c) phase symmetry, and (d) symmetry energy.

3.3.1.

Phase congruency

PC is a dimensionless quantity that measures the agreement of the phase of the Fourier components of the signal (image) being invariant to changes in image brightness or contrast. It differs from gradient-based features because the same relevance is given to all frequency components, independent of gradient magnitude. Hence, an estimate of the noise is removed from the local energy.²⁹

The PC evaluates the local phase information of an image log-Gabor wavelet transform and is computed by

Eq. (24)

$$\mathrm{PC}=\frac{{\sum }_o{\sum }_n{\omega }_o\lfloor {\varphi }_{\mathrm{no}}·{\overline{\varphi }}_o-|{\varphi }_{\mathit{no}}\times {\overline{\varphi }}_o|-T_o\rfloor }{\epsilon +{\sum }_o{\sum }_n{A}_{\mathit{no}}},$$

where $\lfloor .\rfloor$ operator assumes the enclosed quantity as zero when negative, $\epsilon$ is a small positive constant that prevents division by zero, and $A$ is the norm of the phase angle vector $\varphi$.²⁹ $T_o$ is a noise threshold estimated from the response amplitude of the smallest scale filter.²⁹ Moreover, $\varphi$ is given by

Eq. (25)

$${\varphi }_{\mathrm{no}}=[\begin{array}{cc}{M}_{\mathit{no}}^{\text{even}}& M_{\mathit{no}}^{\text{odd}}\end{array}],$$

and the weighted mean phase angle vector ${\overline{\varphi }}_o$ is the unit vector

Eq. (26)

$${\overline{\varphi }}_o=\frac{1}{\sqrt{{\left({\sum }_n{M}_{\mathit{no}}^{\text{even}}\right)}^2+{\left({\sum }_n{M}_{\mathit{no}}^{\text{odd}}\right)}^2}}\left[\sum _n{M}_{\mathit{no}}^{\text{even}}\sum _n{M}_{\mathit{no}}^{\text{odd}}\right],$$

where $M^{\text{even}}$ and $M^{\text{odd}}$ result from the convolution of image $I$ with the even and odd components, respectively, of the log-Gabor kernels. The sigmoid weighting function ${\omega }_o$ is computed by

Eq. (27)

$${\omega }_o(x,y)=\frac{1}{1+\exp \{g[c-s_o(x,y)]\}},$$

where $c$ and $g$ model the sigmoid, and

Eq. (28)

$$s_o(x,y)=\frac{1}{N_n}\left[\frac{{\sum }_n{A}_{\mathit{no}}(x,y)}{\epsilon +A_o^{\max }(x,y)}\right]$$

is the spread of the log-Gabor filters responses.

3.3.2.

Feature type

While higher values of $M^{\text{even}}$ are found at the vessel pixels, $M^{\text{odd}}$ takes higher values at the vessel boundaries and other step/edge locations. To emphasize line-like structures, the even wavelet response is weighted by its odd counterpart. FT was adapted from Ref. 32 and is given by

Eq. (29)

$$\mathrm{FT}=\arctan 2\left[\frac{{\max }_o\lfloor {\sum }_n{M}_{\mathit{no}}^{\text{even}}\rfloor }{{\max }_o({\sum }_n|M_{\mathit{no}}^{\text{odd}}|)}\right],$$

where arctan2 is the four-quadrant arctangent ($\mathrm{FT}\in ]-\pi ,\pi ]$).

3.3.3.

Phase symmetry

PSym is a local contrast-invariant measure of the degree of symmetry of an image.³³ This feature is computed by

Eq. (30)

$$\mathrm{PSym}=\frac{{\sum }_o\lfloor {\sum }_n{M}_{\mathit{no}}^{\text{even}}\rfloor -{\sum }_o{\sum }_n|M_{\mathit{no}}^{\text{odd}}|}{\epsilon +{\sum }_o{\sum }_n{A}_{\mathit{no}}}$$

and was adapted from Refs. 32 and 33.

3.3.4.

Symmetry energy

The last of the local-phase features is the total unnormalized raw symmetry energy³² given by

Eq. (31)

$$\mathrm{SymE}=\sum _o\sum _n[|M_{\mathit{no}}^{\text{even}}|-|M_{\mathit{no}}^{\text{odd}}|-T_o].$$

3.4.

Band-Pass Filter

The PCF image is filtered with a band-pass filter, defined by the log-Gabor radial component as in Eqs. (19) to (21). Several values of central frequency $f_n^0$ were tested. See Table 1 for the chosen parameters. A band-pass filtered PCF can be seen in Fig. 7(a).

Table 1

Parameter values used for feature computation.

Features	Parameters
Expected value fundus	$\rho =7$
Error to local median	$\tau =0.5$
	$w=0.1$
$I^{\text{range}}$, $I^{\mathrm{avg}}$, $I^{\text{entropy}}$, $I^{\mathrm{std}}$	$w=0.03$
${\lambda }^{\max }$, ${\theta }^{\mathrm{var}}$, $\mathrm{\Delta }{I}^{\max }$, ${|\nabla I|}^{\max }$	${\sigma }_n=\{1,4\}$
Phase congruency, feature type, phase symmetry, symmetry energy, $B^{1G}$	$o=\{1,2,\dots ,6\}$
	$n=\{1,2,3,4\}$
Band-pass filter	$n=3$
$B^{\mathrm{avg}}$	$o=\{1,2,\dots ,12\}$
	$m=5$

Equations	Parameters
(11), (33)	$\alpha =6000/(128-1)$
(19)	$\varsigma =0.55$
(20)	$f^{\max }=1/3$
(21)	$f^c=0.45$
	$\upsilon =15$
(27)	$c=0.5$
	$g=10$

Fig. 7

Band-pass (a), average filter-bank (b), and log-Gabor filter bank (c) features computed from the principal component fundus image.

4.

Data

OCT data were gathered from our institution’s database. These OCTs had been acquired resorting to the high-definition spectral-domain Cirrus™ HD-OCT (Carl Zeiss Meditec Inc., Dublin, CA). It allows the acquisition of volumetric data from a region of the human retina of $6000\times 6000\times 2000{\mu \text{m}}^3$, covering a 20 deg field of view of the ocular fundus, with $200\times 200\times 1024$ or $512\times 128\times 1024\text{voxels}$, along the $x$ (nasal-temporal), $y$ (superior-inferior), and $z$ (anterior-posterior) directions, respectively.

Three different datasets—DS1, DS2, and DS3—were used to validate our algorithm.

OCT macular scans of 10 eyes from 10 healthy subjects and 20 eyes from 13 patients diagnosed with type 2 diabetes mellitus [early treatment diabetic retinopathy study (ETDRS) levels 10 to 35] were used for optimization and cross-validation of the classification process. The healthy group contains four volumes of $200\times 200\times 1024\text{voxels}$ and six volumes of $512\times 128\times 1024\text{voxels}$. The diabetic group consists of 11 volumes of $200\times 200\times 1024\text{voxels}$ and 9 volumes of $512\times 128\times 1024\text{voxels}$. This dataset will be referred to as DS1.

Although our objective is the segmentation of the vascular network within the macular region, we also consider the ONH region to allow the comparison with previously proposed methods. With this purpose, a dataset (DS2) composed of ONH-centered OCT scans of 10 eyes from 10 healthy subjects (all acquired with the $200\times 200\times 1024$ protocol) was used.

Finally, dataset DS3 was used to evaluate the robustness of the proposed method when applied to eyes with pathological disorders. It consists of macular OCTs of eight eyes with different pathologic disorders (four OCTs of each protocol).

All 48 PCF reference images were manually segmented pixel-by-pixel by two graders (T.M. and S.S.) who established two ground truths (one per grader) for each image. The first ground truths (T.M.) are used for all SVM training and testing processes (instead of the union or intersection of the two segmentations), following the approach used in Ref. 18. The segmentations of the second grader (S.S.) are used only to assess the intergrader agreement.

5.

Classification

SVM is a supervised-learning algorithm widely used in pattern recognition.^35,36 In this work, pixel-by-pixel classification was performed by SVM. For training, manual segmentations of the vascular network are used to derive the best SVM model.

A C-support vector classification with a radial-basis-function kernel was used,³⁷ where two additional parameters intrinsic to the SVM are required: the parameter $C$ controlling the separability margin of the hyperplane and the parameter $\gamma$ controlling the spread of the kernel. The best $(C,\gamma )$ combination was searched for using a genetic-algorithm aiming for the highest accuracy using cross-validation. The dataset DS1 and the ground truths of the first grader were used for optimization. Because of the large amount of data and the required computing time, we used a twofold cross-validation for the optimization and only a fraction of each image. Specifically, we used 10% of vessel pixels and 10% of nonvessel pixels (both randomly selected) of each image.

All the defined features were used in the SVM. These features were selected from a pool using a forward-selection approach based on the accuracy of the classification. This larger pool of features included all the ones described herewith, as well as variations of these, e.g., different $\rho$, $\tau$, and $w$ [Eqs. (2) and (3)], and additional features like moment invariants-based features^16,38 and semantic/categorical features using spin descriptors,^39,40 to name a few. The segmentations of the first grader and the dataset DS1 were used in this process.

6.

Metrics

Accuracy, specificity, and sensitivity were used for the evaluation of the system performance.

In addition, appropriate metrics were borrowed from Ref. 41. In brief, these are the connectivity C, area A, length L, and CAL. C, A, and L are given by

Cohen’s $\kappa$ (Ref. 42) was computed as a metric of segmentation agreement. It is defined by

7.

Results

Each of the features used by the supervised-classification process requires the specification of working parameters. These parameters were defined based on the system characteristics (e.g., scanning protocol), the forward-selection method (features), or values previously defined in the literature. The parameters used in this work can be found in Table 1.

It is important to note that no postprocessing was applied. All automatic segmentations are the direct result of the classification process.

The classification was validated using a 10-fold cross-validation approach on the dataset DS1 with manual segmentations from the first grader. Each of the 10 groups was composed of macular scans from one healthy and two diabetic subjects. OCTs acquired with different protocols were also distributed equitably throughout the groups. Quantitative results are summarized in Table 2. For visual inspection, the cases with worst and best performance (for the accuracy and CAL metrics) are shown in Fig. 8.

Table 2

Results from the 10-fold cross-validation on dataset DS1, healthy and diabetic subjects. The minimum (MIN), maximum (MAX), average (AVG), and standard deviation (SD) values for each segmentation performance metric.

Fig. 8

Principal component fundus (top row), and manual (middle row) and automatic (bottom row) segmentations from optical coherence tomography volumes with the highest (a) and lowest accuracy (b), and the highest (c) and lowest CAL (d).

The achieved values demonstrate the viability of the proposed method. The high specificity (0.994) shows that there is little oversegmentation, percentagewise. On the other hand, undersegmentation is a more prominent problem, as demonstrated by the lower sensitivity (0.825). By visual assessment of Fig. 8, one can say that the undersegmentation is associated with the misclassification of pixels on low-caliber vessels. Considering the application of subsequent algorithms (e.g., 3-D vascular segmentation and computation of vascular network descriptors), a specificity close to 1 (associated with a high sensitivity) is particularly important, as the postprocessing stage would be less demanding.

The differences between the healthy and diabetic groups are small. However, the algorithm performed better on the healthy group.

A new training was performed using all OCT scans in DS1. The model was then used to segment the OCT scans in DS2, the ONH OCTs, and in DS3, the pathological cases. For the ONH dataset, the optic disk was manually segmented and discarded in metric computation. Quantitative results are shown in Tables 3 and 4. For qualitative assessment, Fig. 9 shows a few examples of vascular segmentation on the ONH OCTs, and Fig. 10 shows all segmentation results of the pathological cases. Although the segmentation was developed mainly for OCTs of healthy eyes and eyes close to the healthy condition, the results on the pathological cases serve to demonstrate the robustness of the proposed method in rather extreme conditions.

Table 3

Test results for the optic nerve head dataset (N=10). The minimum (MIN), maximum (MAX), average (AVG), and standard deviation (SD) values for each segmentation performance metric.

Table 4

Test results for the pathology dataset (N=8). The minimum (MIN), maximum (MAX), average (AVG), and standard deviation (SD) values for each segmentation performance metric.

Fig. 9

Examples of optic nerve head optical coherence tomography fundus images and the automatic segmentation results. The gray region (from the manual segmentation of the optic nerve head) was not considered.

Fig. 10

Optical coherence tomography fundus images and automatic segmentation results for cases with pathological disorders: (a) full-thickness macular hole, (b) branch retinal arterial occlusion, (c) cystoid macular edema, (d) and (e) age-related macular degeneration, (f) vascular tortuosity, (g) proliferative diabetic retinopathy, and (h) diabetic macular edema.

The proposed approach was able to outperform the one published in Ref. 26 for segmentation on the ONH region, where a specificity of 88% and a sensitivity of 85% were achieved. Although the dataset is not the same, it has the same fundamental characteristics: the OCT scans were acquired with Cirrus HD-OCT using the same acquisition protocol ($200\times 200\times 1024$), all OCTs were centered on the ONH region, the number of cases is similar ($N=11$), and all eyes are from healthy subjects. Furthermore, in our work, the optimizations and training processes were performed on OCT scans from the macular region. Therefore, the process could benefit from parameter optimization and training on ONH-centered scans.

In the pathological cases, the specificity is lower, which translates to higher oversegmentation. However, while no training was performed on the pathological cases other than the diabetic ones, the results achieved on several pathologies are similar to the results achieved with the healthy and the diabetic subjects (accuracy $>97\%$ and $\kappa >0.81$), which attests to the robustness of the proposed method.

The agreement between the two graders was also evaluated. The manual segmentations of the second grader (S.S.) were tested against those of the first grader (T.M.). Results can be found in Table 5. These confirm the good performance of the algorithm. For the majority of the metrics, the average values for the automatic segmentations on DS1 and DS2 are greater than or equal to those for the second grader. Contrarily, as expected, for DS3, the metrics (although similar) are lower than those for the comparison between graders.

Table 5

Intergrader agreement (N=48). The minimum (MIN), maximum (MAX), average (AVG), and standard deviation (SD) values for each segmentation performance metric.

The total execution time for the segmentation process (OCT fundus reference computation, features computation, and SVM classification), using a nonoptimized MATLAB^® (The MathWorks Inc., Natick, MA) implementation, was $65.2\pm 1.2\mathrm{s}$ ($N=15$) and $108.2\pm 3.8\mathrm{s}$ ($N=15$) for the $200\times 200\times 1024$ and $512\times 128\times 1024$ protocols, respectively. The system hardware used was an Intel^® Core™ i7-3770 CPU (Intel Corporation, Santa Clara, California) at 3.4 GHz.

8.

Discussion and Conclusions

To our knowledge, the work presented herein is the first attempt to automatically segment the vascular network of the macular region. Our findings showed that it is able to work with standard OCTs of both healthy and diseased retinas. The proposed method achieved good results for both the macular and ONH regions even with no training with ONH OCTs.

The proposed models (trained only on healthy and diabetic subjects) already show robustness on the pathological cases, as the results do not differ substantially from DS1 results. Nevertheless, with proper training on pathological data and given enough pathological cases, the SVM would be able to create even more robust models, which would improve the classification performance on these cases.

Improvements are still possible and the 2-D segmentations could benefit from a postprocessing stage. The spectrum of features presented here can be useful for region growing based on the automatic segmentations and subsequent deletion of small incorrectly segmented regions. For the pathological cases, the classification specificity is lower, which leads to a more complex postprocessing stage than for DS1 and DS2. Additional research will be performed on these issues.

With the proposed method, additional algorithms of coregistration⁴³ and computation of vascular network descriptors,^44,45 already described for other imaging modalities, can now be applied to OCT. For most multimodal coregistration algorithms, even for the pathological cases, we are particularly convinced that the proposed segmentation is a valuable tool. Therefore, multimodal imaging and studies of disease progression (even in extreme cases) are a major application area of the proposed method.

The representation of blood vessels on OCT fundus reference images depends on the presence of hemoglobin. Therefore, the method is not able to segment fully occluded vessels. However, for reperfused or partially occluded vessels, segmentation is still possible [Fig. 10(b)]. This led to yet another area of application now being pursued by our research group: the discrimination between perfused and occluded vessels from OCT data.⁴⁶

Finally, this method establishes a good starting point toward the fully automatic 3-D segmentation. The first tests on 3-D segmentation have been conducted in Ref. 25.