Abstract

Spectral-Domain Optical Coherence Tomography (SD-OCT) is a widely used interferometric diagnostic technique in ophthalmology that provides novel in vivo information of depth-resolved inner and outer retinal structures. This imaging modality can assist clinicians in monitoring the progression of Age-related Macular Degeneration (AMD) by providing high-resolution visualization of drusen. Quantitative tools for assessing drusen volume that are indicative of AMD progression may lead to appropriate metrics for selecting treatment protocols. To address this need, a fully automated algorithm was developed to segment drusen area and volume from SD-OCT images. The proposed algorithm consists of three parts: (1) preprocessing, which includes creating binary mask and removing possible highly reflective posterior hyaloid that is used in accurate detection of inner segment/outer segment (IS/OS) junction layer and Bruch’s membrane (BM) retinal layers; (2) coarse segmentation, in which 3D curvelet transform and graph theory are employed to get the possible candidate drusenoid regions; (3) fine segmentation, in which morphological operators are used to remove falsely extracted elongated structures and get the refined segmentation results. The proposed method was evaluated in 20 publically available volumetric scans acquired by using Bioptigen spectral-domain ophthalmic imaging system. The average true positive and false positive volume fractions (TPVF and FPVF) for the segmentation of drusenoid regions were found to be 89.15% ± 3.76 and 0.17% ± .18%, respectively.

1. Introduction

Age-related Macular Degeneration (AMD) is the most common degenerative eye disease that affects the central portion of the retina, known as the macula, and is the leading cause of vision loss and blindness in patients who are aged 65 and over in the developing countries [1]. There are two types of AMD: dry (atrophic or nonneovascular) and wet (neovascular or exudative) [2]. Dry form of AMD affects approximately 80–90% of individuals with AMD [3] and in this type of macular degeneration, extracellular small white or yellowish deposits (made up of lipids, a type of fatty protein), called drusen, accumulated between the Retinal Pigment Epithelium (RPE) and the inner collagenous layer of Bruch’s membrane (vitreous lamina). Drusen are produced when the photoreceptors of the eye drop off older parts of the cell and could cause a protrusion of RPE over Bruch’s membrane which leads to deterioration or degeneration of the macula over time [4].

Drusen in the macula are a common early sign of AMD that can be defined as hard and soft. Hard drusen are small and distinct lesions with sharp borders (this type of drusen may not cause vision problems for a long time and is normal with aging, and most people over 40 have some hard drusen). In contrast, soft drusen are those with large lesions (clusters close together) with indistinct borders that have been recognized as intermediate to advanced form of AMD [5].

Most patients with macular degeneration have the dry form of the disease. However, the dry form of macular degeneration can lead to the more damaging wet form that usually leads to more serious vision loss [2]. The Age-Related Eye Disease Study (AREDS) reported that certain antioxidant vitamins and zinc delay the progression of advanced dry AMD or vision loss by about 25 percent over a six-year period [6]. The early diagnosis of the patients at higher risk for advanced AMD allows earlier protective intervention and preventive treatment to reduce severe vision loss. Thus, it is more important to identify the main changes in drusen and the RPE morphology to inspect the progression of nonneovascular AMD [7]. As shown in Table 1, AREDS defines four categories depending on the stages of AMD [8].

The previous studies on the color fundus (CF) images have shown that there is strong correlation between the maximum drusen size and total drusen area, with the risk of progression of AMD to its advanced form [1, 8]. Many different approaches and algorithms for automatic soft and hard drusen grading and quantification from CF images have been developed [912]. However, detecting and locating the drusen in a color retinal image are difficult to measure because of their differences in size and shape and also irregular and blurred boundaries.

Spectral-Domain Optical Coherence Tomography (SD-OCT) is a widely used interferometric diagnostic technique in ophthalmology that provides in vivo novel information of depth-resolved inner and outer retinal structures. The high resolution, sensitivity, and speed of SD-OCT provides a valuable opportunity for detection and volume measurement of drusen from SD-OCT images compared to color fundus photographs [13] and can be used as a method of following the AMD development and monitoring retinal changes over the time [8, 14, 15]. However, manual segmentation of macular drusen area and volume is time consuming and labor intensive, which limits its use in large-scale population. So, the automatic segmentation of drusen becomes more appealing to process SD-OCT data more efficiently.

OCT can provide useful information about drusen and results in more accurate differentiation of drusenoid deposits, including cuticular drusen, soft drusen, and subretinal drusenoid deposits [16]. Soft drusen are yellow-white dome shaped elevations that their central portion appear lighter than its edge and are typically 63 to ≥1000 μm in diameter. Smaller soft drusen have a little effect on the RPE and IS/OS junction layer morphology. Cuticular drusen appear as a round and punctate, with saw tooth pattern in an SD-OCT and are typically 50 to 75 μm in diameter [16]. But, OCT subretinal drusenoid deposits are seen in the subretinal space, in similar size to soft drusen, above the RPE. Soft and cuticular drusen, in contrast, are under the RPE. They often have a more punctate appearance and are interconnected so that their sizes range from 25 to >1000 μm [16].

As shown in Figure 1, small and medium-size drusen may cause subtle elevation of discreet areas of RPE with variable reflectivity, indicating the variable accumulations of the underlying material. Greater elevation of RPE may be seen as the result of larger confluent drusen or Drusenoid Pigment Epithelium Detachment (DPED), with a hyporeflective or medium-reflective deposit separating the RPE from the underlying Bruch membrane (BM) [17].

Although a number of fully/semiautomatic techniques have been proposed to detect drusen from OCT images (Table 2), only the accuracy of a few methods has been validated and they mainly only identify the absence or presence of drusen. The challenging issues associated with drusen detection from OCT images include the variation of texture and intensity distribution (hyporeflective or medium-reflective) of drusen which prevents accurate segmentation of the RPE and BM for quantification of drusen complex (the distance between the abnormal RPE and the normal estimated RPE floor), and methods which are based on accurate detection of RPE layer [18, 19] may fail due to difficulty in RPE layer segmentation even in normal patients, because the RPE layers are often contiguous with the inner segment/outer segment (IS/OS) retinal layer.

In this paper, we tackle the above-mentioned challenges and present a novel fully automated drusen segmentation method in SD-OCT images. For this purpose, in order to prevent inaccurate detection of retinal layers, we preprocess OCT images and remove highly reflective posterior hyaloid or other reflective layers beneath the BM layer, then we use curvelet transform and modify curvelet coefficients to enhance and fill the possible appeared vertically black shadowing effects due to blood vessels. These possible gaps due to large vessels may degrade graph-based segmentation of IS/OS and BM layers. The areas located between the shifted extracted IS/OS junction layer and BM pixels by the use of graph theory from modified image are marked as drusen.

The paper is organized as follows. Section 2 provides an introduction to 3D digital curvelet transform (3DCUT). In Section 3 our proposed method is described and the results and performance evaluation are presented in Section 4. Finally, the conclusions are given in Section 5.

2. 3D Digital Curvelet Transform

The curvelet transform is a time-frequency analysis of images that gives a sparse representation of objects. The basis elements of this transform have good directional selectivity and are highly anisotropic [20]. The directional selectivity of curvelets and spatially localized property of each curvelet can be utilized to preserve the image features along certain directions in each subband. Following this reasoning, curvelets are appropriate basis elements (atoms) for representing d-D objects which are smooth apart from singularities along smooth manifolds of codimension 1 [21].

Although the direct analyzing of 3D data as a volume and considering the 3D geometrical nature of the data are computationally expensive, it has been shown that 3D analysis of 3D data outperforms 2D slice-by-slice analyzing of data [22]. The parabolic scaling, good directional selectivity, tightness and sparse representation properties of 3D curvelet singularities, provide new opportunities to analyze and study large data-sets in medical image processing. 3D curvelet elements are plate-like shapes of in two directions and width about in the orthonormal direction which are smooth within the plate and oscillate along the normal direction of the plate.

The 3D discrete curvelet transform is an extension to the 2D curvelet transform proposed in [21]. In 2D, the curvelet family is generated by translation and rotation of the basic element :where is rotation matrix with angle and is given by its Fourier transform . The -scaled window is a polar wedge in frequency domain that is used for building curvelet functions [23]. Sinceusing a suitable sampling at the range of scales , orientations ( denote the smallest integer being greater/smaller than or equal to ), and locations , the curvelet coefficients are defined as

To have Cartesian arrays instead of the polar tiling of the frequency plane, the Cartesian window is defined [23] that determines the frequencies in the trapezoid . So, the Cartesian counterpart of the coefficients can be defined bywhere , and and .

By defining regular rectangular grid instead of tilted grid, Cartesian curvelets are defined as

of formula (4) has been replaced by that is being taken on values on a rectangular grid.

In 3D, similar to 2D case, the Cartesian window is defined by that isolates the frequencies in the truncated pyramidWith the angles and the 3D shear matrix is defined as , whereIn 3D, every Cartesian corona has six components, one for each face of the unit cube. Each component is regularly partitioned into wedges with same volume as shown in Figure 2.

The curvelet functions in the cone are given bySo, its Fourier transform would bewhere Analogously as in (5), the curvelet coefficients are given byIn this paper we use a new implementation of the 3D Fast Curvelet Transform (3DFCT) [24, 25] with strong directional selectivity property at the finest scale that has a reduced redundancy factor than the wrapping-based implementation proposed in CurveLab Toolbox [20, 21]. According to this implementation, 3D curvelet coefficients are obtained as follows:(1)Cartesian coronization is performed that decomposes the object into dyadic coronae based on concentric cubes. Each corona is subdivided into trapezoidal regions conforming the usual parabolic scaling as shown in Figure 2.(2)The 3D coefficients are obtained by applying an inverse 3D FFT to each wrapped wedge as shown in Figure 2, which appropriately fits into a 3D rectangular parallelepipeds of dimensions ~ centered at the origin.

3. Proposed Method

3.1. Preprocessing

The posterior hyaloid membrane separates the vitreous from the retina. A Posterior Vitreous Detachment (PVD) is a condition of the eye, which is common in over 75% of patients who are age 65 and older, in which the vitreous membrane is separated from the retina [26]. As shown in Figure 3, the presence of detached posterior hyaloid (thin layer above retina) in some B-scans of 3D OCT volumetric data, creates an abnormal layer that may be seen as a thin hyperreflective layer above the internal limiting membrane (ILM). Some previous algorithms [27, 28] that assume ILM to be the first hyperreflective layer may determine this abnormal detached vitreous membrane as the retinal boundary, which results in an inaccurate segmentation of the remaining layers, also overestimation of retinal thickness.

To overcome this problem, each OCT image is thresholded with an optimal threshold selected efficiently for each image based on entropy-based thresholding algorithm [29], which takes into account the spatial distribution of gray levels in order to coarsely extract the Region of Interest (ROI) from the nonretinal layer background.

Some possible connected pixels due to highly reflective posterior hyaloid in thresholded image can be removed by applying length filtering and removing elongated short components. In the extracted ROI, possible holes near ILM are also removed by dilating and eroding (morphological operators) with the disk shaped structuring element of radius 5.

In some cases upper band of ROI is segmented with some kind of bumps (Figure 4(b), right). In order to segment this layer smoothly, we fit a third order polynomial to the upper pixels of ROI region. We use this boundary to limit the search region in segmentation process of IS/OS and BM and differentiate them from ILM (3 distinct hyperreflective layers with similar characteristics). Once the binary mask is generated, 10 pixels below the upper detected boundary in created binary image is used to separate inner and outer portion of retina and reduce the search range which results in an accurate and fine segmentation of BM and IS/OS junction layer by the use of graph theory.

Figure 4(d) shows the created binary mask in two sample B-scans contain posterior hyaloid and bumps in detected upper boundary of ROI.

3.2. Drusenoid Region Segmentation

The appeared vertically black shadowing effects due to blood vessels may degrade the continuity of IS/OS junction layer, and methods based on thresholding and graph theory may fail for accurate detection of these layers in presence of large vessels. So, in order to enhance IS/OS layer and fill these gaps we modify curvelet coefficients in OCT images. Since the curvelet coefficients have a sparse distribution and the bright objects get large coefficients in the curvelet domain in order to enhance the bright hyperreflective regions in OCT images, similar to one defined by Starck et al. [30], we take curvelet transform of OCT image and modify DCUT coefficients using function :In this equation , where is the maximum curvelet coefficient of the relative band. Then, we reconstruct the enhanced image from the modified curvelet coefficients by applying inverse DCUT. Figure 5(b) shows the reconstructed image from modified curvelet coefficients.

To delineate the BM, we use modified reconstructed image with upper masked region using the extracted binary mask in previous section (Figure 5(c)). By using this image that the upper cropped region is set to zero, the BM layer would be the first light-to-dark boundary (vertically). Then, we apply the graph theory based segmentation method in [31] with an adjacency matrix that its weights are calculated using the following equation:In this equation, is the weight assigned to the edge connecting nodes and , and are the vertical dark-to-light gradients of the image at nodes and , and are the vertical light-to-dark gradients of the image at nodes and and is the Euclidian distance from node to node , and 1 × 10−5 is the minimum weight of the graph.

The norm notation indicates the normalization of the value to range from to .

For IS/OS junction boundary, we use modified reconstructed image after setting the upper cropped regions to maximum intensity of image (Figure 5(d)). By using this image, the IS/OS layer is the first vertical dark-to-light boundary. For this image, the adjacency matrix is constructed by calculating the weights as follows:Figures 6(b) and 6(c) illustrate the results of applying the mentioned graph-based method on the masked images for a sample B-scan (Figure 6(a)).

In order to estimate RPEDC (RPE + drusen complex) thickness, at first we apply two 2nd order polynomial fitting functions to the detected IS/OS and BM (Figure 6(d)), then the mean difference of these two 2nd order polynomial curves is calculated. The calculated mean value is used to estimate RPEDC thickness (Figure 6(e)) and to shift upper detected IS/OS pixels toward detected BM (Figure 6(f)).

For extraction of drusenoid regions, we produce an image that in each column of this image if row location (in MATLAB coordinate system) of IS/OS pixel + mean difference is less than or equal to the row location of detected BM we set the row location of shifted IS/OS (row location of IS/OS pixel + mean difference) and row location of BM in that column to 1; otherwise this column is set to zero (Figure 6(f)). The created holes between connected layers (shifted IS/OS and BM) are filled by performing a flood-fill morphological operation on the pixels of the binary image (Figure 6(g)). Then we remove miss extracted liny shape structures by morphological erosion and dilation of image with disk shape structure element (Figure 6(h)). Figure 6(i) shows the final results of our proposed method.

A sketch of the proposed algorithm after detecting IS/OS and BM layers with the use of graph theory is as in Algorithm 1.

Input: I% I is an image , with extracted locations of IS/OS and BM (Figure 6(d))
Output: DR% Create the image DR, for extraction of drusenoid regions
DR = zeros(size(I));
For
;
= Find ;
If + mean_diff-max(i)<=0% mean_diff is obtained from I (Figure 6(e))
DR(min(i) + mean_diff,l) = 1;
DR = 1;
End
End

Figure 7 illustrates the general block diagram of our proposed approach.

4. Results

For this study we use 20 volumetric scans acquired by using Bioptigen spectral-domain ophthalmic imaging system (the dataset is publicly available at http://people.duke.edu/~sf59/Chiu_IOVS_2011_dataset.htm) [31]. Two independent expert graders performed manual segmentation of the drusenoid regions using the 3D Slicer software (the software is freely available at https://www.slicer.org/). We then performed automatic segmentation using the algorithm described earlier, which was implemented in MATLAB.

To test the accuracy of proposed algorithm, automatic segmentation results for 220 B-scans across 20 patients were compared against the available results from two graders.

Segmentation results are quantitatively measured using Dice Coefficient (DC) and true positive and false positive volume fractions (TPVF and FPVF) [32] are defined as follows:where and are the extracted volume in the automatic segmented image and in the reference ground truth image. The maximum possible value for DC is 1 (indicating a perfect match between the result of the algorithm and ground truth) and DC is zero for a slice when no drusenoid region is detected or when there is no overlap between and .

The mean and standard deviation of DC for all the volumes with ground truth taken as Graders 1 and 2 and union of Graders 1 and 2 are mentioned in Table 3.

TPVF also indicates the fraction of the total amount of tissue inside the reference delineation, and FPVF denotes the amount of tissue falsely identified. They are defined as follows:where denotes the volume in the whole image between ILM and BM.

In order to show the effect of applying 3D curvelet transform in the proposed method in this paper, we can implement the different stages of proposed segmentation algorithm without applying 3D curvelet transform and compare the results. As an example, Figure 8 shows the effect of applying 3D curvelet transform in accurate detection of BM. In comparison to the proposed method in [18], which is based on assuming the RPE layer to be constant in thickness (20 μm) and using the threshold-based method to identify RPE layer, our proposed method can robustly identify RPE and IS/OS layer for drusen quantification. Figure 9 illustrates the final results of proposed drusenoid region segmentation method in this paper for different B-scans.

5. Conclusion

SD-OCT provides high axial resolution and dense 3D imaging of the pathological changes in the eye, which makes it an appropriate device for detection of abnormal elevations of the RPE and IS/OS layer in various pathologies such as drusenoid pigment epithelial detachment. Previous methods that are based on difference between detected normal and abnormal RPE layer may fail due to difficulty in RPE layer segmentation even in normal patients because the RPE layer are often contiguous with the IS/OS retinal layer.

This paper presents a new fully automatic method to segment and quantify the drusen volume from OCT images, in order to provide a metric that can assist clinicians in monitoring the progression of AMD. The previous studies show that there are good agreement of drusen diameter and mean drusen area on SD-OCT with those identified on color fundus photographs. In this base, drusen detection using 3D SD-OCT may be a potentially useful alternative method to drusen assessment by human graders using color fundus photographs.

In order to prevent inaccurate detection of retinal layers by graph-based intraretinal layer segmentation methods, we preprocessed OCT images and removed highly reflective posterior hyaloid (or other reflective layers above RPE), and then we applied 3D curvelet transform and modified curvelet coefficients to enhance and fill the possible appeared vertically black shadowing effects due to blood vessels. These possible gaps due to large vessels may degrade graph-based segmentation of IS/OS and BM layers.

The shifted upper detected IS/OS pixels that are above BM are considered as a possible candidate region and dome shaped regions are considered as a drusenoid regions. The proposed method for drusen segmentation is not able to identify small size drusen (approximately less than 8 pixels). The information about the drusenoid area can be better assessed by creating OCT fundus images as an enface projection image from the SD-OCT dataset and registering these images to corresponding color fundus photographs and mapping the delineated drusen area on each color fundus photograph.

Competing Interests

The authors declare that they have no competing interests.