D Curvelet-Based Segmentation and Quantification of Drusen in Optical Coherence Tomography Images

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 inmonitoring the progression of Age-relatedMacularDegeneration (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 algorithmwas 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 proposedmethod 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.


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 Many medium-sized drusen or one or more large drusen in one or both eyes Degeneration of the deepest cells of the retina (in one eye only) or abnormal and fragile blood vessels known as choroidal neovascularization under the retina (wet form) 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 [9][10][11][12].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 mediumreflective) 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 Gabor filter and wavelet analysis [33] Statistical structural information-based texture analysis [34] The Artificial Neural Network (ANN) based supervised classifiers [35] Expectation Maximization (EM) based unsupervised classifiers [36] Thresholding-based segmentation Probability based thresholding [37] Otsu method for background selection and intensity thresholding [38] Histogram based adaptive local thresholding [39] Edge based thresholding [11] Clustering-based segmentation Gaussian template matching [43] Gradient based segmentation [44] OCT image analysis techniques

Graph-based segmentation Graph and dynamic programming [31]
A graph-based multilayer segmentation approach [45] Active contour-based segmentation Anisotropic noise suppression and deformable splines [46] Based on active contours and Markov random fields [47] Polynomial curve-fitting based technique Local convexity condition and fitting second-or fourth-order polynomials [48,49] Based on the distance between the abnormal RPE and the normal RPE floor [19,50] 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.

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 datasets in medical image processing.3D curvelet elements are plate-like shapes of 2 −/2 in two directions and width about 2 − in the orthonormal direction which are smooth within the plate and oscillate along the normal direction of the plate.
In 3D, similar to 2D case, the Cartesian window Ũ () is defined by Ũ () = Ũ ( 1 ,  2 ,  3 ) that isolates the frequencies in the truncated pyramid With the angles  , and  , the 3D shear matrix is defined as ), where In 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.
So, its Fourier transform would be φ,,, = exp (− ⟨ b,, where φ,0,0,0 () = Ũ ().Analogously as in (5) In 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 ∼ (2  , 2 /2 , 2 /2 ) centered at the origin.

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.

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  = 0.1, 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] Posterior hyaloid  with an adjacency matrix that its weights are calculated using the following equation: + norm (− ( DL  +  DL  ) , 0, 0.5) In this equation,   is the weight assigned to the edge connecting nodes  and b,  DL  and  DL  are the vertical darkto-light gradients of the image at nodes  and b,  LD  and  LD  are the vertical light-to-dark gradients of the image at nodes  and  and   is the Euclidian distance from node  to node , and  min = 1 × 10 −5 is the minimum weight of the graph.The norm notation norm(, , ) 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: + norm (−  , 0, 0.5) +  min 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)).
Modify curvelet coefficients with 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.
Figure 7 illustrates the general block diagram of our proposed approach.

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]   and DC is zero for a slice when no drusenoid region is detected or when there is no overlap between  seg and  ref .
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  .V 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.

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.

Spatial histogram and similarity based classification [ 40 ]
Distance based clustering [41] Wavelet based feature extraction & SVM based classification [42] Edge and template matching

Figure 4 :Figure 5 :
Figure 4: Created binary mask to accurate segmentation of BM and RPE.First and second columns, respectively, show the results for a sample B-scan with posterior hyaloid and a sample data with bumps in detected ROI region.(a) Original image.(b) Thresholded + length filtered + filled ROI image.(c) Segmented upper band pixels of ROI.(d) Created binary mask by fitting 3rd-order polynomial to upper band pixels of ROI.

Figure 6 :
Figure 6: Drusenoid region segmentation results.(a) Original noisy image, (b) detected IS/OS junction boundary, (c) detected BM layer, (d) extracted BM and IS/OS boundaries, (e) 2nd-order polynomial fitted lines to extracted boundaries in order to calculate the mean difference of two boundaries, (f) created image by shifting the upper detected IS/OS junction pixels, (g) filled image, (h) reconstructed image by applying morphological erosion and then dilation operators, and (i) segmented drusenoid region.
OCT image (I)Apply lengthfiltering  and remove elongated short components Dilate and erode D with the disk shaped structuring element to fill possible holes Take curvelet transform of image, and collect curvelet coefficients Take inverse curvelet transform of modified coefficients Use upper detected pixels of ROI and limit the search region for graph-based segmentation of IS/OS junction layer and BM layer from reconstructed image Apply 2nd-order order polynomial fitting functions to the detected IS/OS and BM Calculate the mean difference (m) of these two 2nd-order polynomial curves Shift upper detected IS/OS by m pixels toward detected BM Fill the holes in produced image by shifting IS/OS toward BM Remove miss extracted line shape structures by morphological erosion and dilation of image with disk shape structure element Apply threshold T is selected threshold based on entropy-based thresholding

FitFigure 7 :
Figure 7: The block diagram of the proposed method.

Figure 8 :
Figure 8: The results of detected BM by applying the graph-based segmentation method on (a) modified image by curvelet transform and (b) original image without modification.

Figure 9 :
Figure 9: Experimental results for three examples of drusen segmentation.The first column shows the original image and the 2nd column shows the final fine segmentation results.

Table 2 :
Available image analysis techniques for AMD.
I% I is an image  × , with extracted locations of IS/OS and BM(Figure 6(d)) where  seg and  ref 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) Input:

Table 3 :
Comparison of automatic versus manual results of drusenoid region segmentation.