Issuu

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication. The final version of record is available at http://dx.doi.org/10.1109/TBME.2017.2734058 JOURNAL OF LATEX CLASS FILES, VOL. XX, NO. X, XXXX

Fully-Automated Segmentation of Fluid/Cyst Regions in Optical Coherence Tomography Images with Diabetic Macular Edema using Neutrosophic Sets and Graph Algorithms Abdolreza Rashno, Student Member, IEEE, Dara D. Koozekanani, Member, IEEE, Paul M. Drayna, Behzad Nazari, Saeed Sadri, Hossein Rabbani, Senior Member, IEEE, Keshab K. Parhi, Fellow, IEEE,

Abstract—This paper presents a fully-automated algorithm to segment fluid-associated (fluid-filled) and cyst regions in optical coherence tomography (OCT) retina images of subjects with diabetic macular edema (DME). The OCT image is segmented using a novel neutrosophic transformation and a graph-based shortest path method. In neutrosophic domain, an image g is transformed into three sets: T (true), I (indeterminate) that represents noise, and F (false). This paper makes four key contributions. First, a new method is introduced to compute the indeterminacy set I, and a new λ-correction operation is introduced to compute the set T in neutrosophic domain. Second, a graph shortest-path method is applied in neutrosophic domain to segment the inner limiting membrane (ILM) and the retinal pigment epithelium (RPE) as regions of interest (ROI) and outer plexiform layer (OPL) and inner segment myeloid (ISM) as middle layers using a novel definition of the edge weights. Third, a new cost function for cluster-based fluid/cyst segmentation in ROI is presented which also includes a novel approach to estimating the number of clusters in an automated manner. Fourth, the final fluid regions are achieved by ignoring very small regions and the regions between middle layers. The proposed method is evaluated using two publicly available datasets: Duke, Optima, and a third local dataset from the UMN clinic which is available online. The proposed algorithm outperforms the previously proposed Duke algorithm by 8% with respect to the dice coefficient and by 5% with respect to precision on the Duke dataset, while achieving about the same sensitivity. Also, the proposed algorithm outperforms a prior method for Optima dataset by 6%, 22% and 23% with respect to the dice coefficient, sensitivity and precision, respectively. Finally, the proposed algorithm also achieves sensitivity of 67.3%, 88.8% and 76.7%, for the Duke, Optima, and the UMN datasets, respectively. Index Terms—Fluid/cyst segmentation, graph theory, neutrosophic set, optical coherence tomography, diabetic macular edema.

A. Rashno is a student at Isfahan University of Technology, Isfahan 84156, Iran and visiting research student at University of Minnesota, Minneapolis, MN 55455, USA (e-mails:ar.rashno@gmail.com). B. Nazari and S. Sadri are with the Department of Electrical and Computer Engineering, Isfahan University of Technology, Isfahan 84156, Iran (emails:ar.rashno@gmail.com, nazari@cc.iut.ac.ir and sadri@cc.iut.ac.ir). D. D. Koozekanani and P. M. Drayna are with the Department of Ophthalmology and Visual Neurosciences, University of Minnesota, Minneapolis, MN 55455 USA (emails: dkoozeka@umn.edu and drayn006@umn.edu). H. Rabbani is with the Department of Biomedical Engineering, Medical Image and Signal Processing Research Center, Isfahan University of Medical Sciences, Isfahan, Iran (e-mail: rabbani.h@ieee.org). K. K. Parhi is with the Department of Electrical and Computer Engineering, University of Minnesota, Minneapolis, MN 55455, USA (email: parhi@umn.edu).

I. I NTRODUCTION PTICAL coherence tomography (OCT) is a non-invasive and non-contact imaging method with extensive clinical use in ophthalmology. It uses optical technology to create tomographic images with variable scan rates and resolutions, and is used to create cross-sectional images of ocular tissues, including the retina [1], [2]. It is extensively used clinically for the diagnosis and follow-up of patients with DME and age-related macular degeneration (AMD). OCT images allow detection and quantitative assessment of retinal abnormalities [3]. The macula is the central part of the retina and is critical for good vision. Diabetic macular edema (DME), manifested by fluid cysts within the retina and retinal thickening, is caused by fluid leakage from damaged macular blood vessels. This is the most common cause of vision loss among workingaged adults in the United States. OCT images allow very sensitive detection and quantitative assessment of these fluid cysts and retinal thickening [4], [5]. Ophthalmologists evaluate the severity of DME using retinal thickness maps to indirectly measure the intra-retinal fluid/cyst regions in OCT images. Manual segmentation of fluid/cyst regions is a time consuming task and is prone to human errors [6]. Currently, very few automated methods for fluid/cyst segmentation exist; some of these are reviewed in Section II. Neutrosophy is a branch of philosophy which studies the nature and scope of the neutralities and their interactions which is the basis of neutrosophic logic and neutrosophic (NS) set [7]. This theory was applied for image processing first by Guo et al. [8] and it has subsequently been successfully used for other image processing operations including image segmentation [8]–[11], image thresholding [12], medical image segmentation [13] and edge detection [14]. Also, NS has been adapted for data and image clustering as well [15]. Graph theory is one of the powerful tools for image segmentation due to the benefits of mapping the image pixels (voxels) and relationship between them onto a graph. Graph theory based image segmentation makes use of techniques such as minimal spanning tree, graph cut with cost function, graph cut on Markov random field models, shortest path methods and random walker methods [16]. Graph-based image segmentation methods have been applied in OCT segmentation; some

prior works include: automated layer segmentation [2], [17]â€“ [19], optic disc segmentation [20] and drusen segmentation [21]. The main contributions of this paper that lead to a fullyautomated method for the segmentation of fluid-associated and cyst regions in two-dimensional (2D) OCT of DME subjects are summarized as follows: 1) The OCT images are first transformed into three sets T (true), I(indeterminacy) and F (false) in the NS domain and denoised simultaneously. This is done by the novel definition of indeterminacy set and proposed Îťcorrection operation which is robust to the structures and intrinsic properties of the OCT images. The traditional NS domain processing, developed for natural images, is not directly applicable to OCT segmentation due to its layered structure. The proposed approach is the first approach to make use of NS logic for OCT image segmentation. Set T is the output of this step and is used for all subsequent steps. 2) The outer layers inner limiting membrane (ILM) and retinal pigment epithelium (RPE) are segmented as region of interest (ROI) and the middle layers outer plexiform layer (OPL) and inner segment myeloid (ISM) are segmented using the proposed graph shortest path based method in the NS domain. With respect to graph construction, new definitions of edge weights are introduced for layer segmentation. Note that middle-layer segmentation is a key step of the algorithm to reduce the number of false positives which is used in the postprocessing step. 3) A new method in the neutrosophic domain is proposed for fluid/cyst segmentation. Then, the appropriate number of clusters is computed automatically from the estimated fluid/cyst regions and is used by the proposed clustering scheme for segmentation. In the clustering method, a new cost function is defined and minimized to obtain the segmented image. 4) In the post-processing step, both the segmented regions between middle layers and small segmented regions are ignored. The rest of this paper is organized as follows: Section II describes related works. The neutrosophic set approach is reviewed in Section III. Proposed algorithms are presented in Section IV. Experimental setup and results are described in Sections V and VI, respectively. Statistical comparison and discussion are presented in Section VII. Finally, conclusions and future directions are discussed in Section VIII. A preliminary version of this work has been reported in [22]. II.

RELATED WORKS

For segmentation and quantification of fluid/cyst features in OCT images, several methods have been proposed. A semiautomated method for fluid-filled region segmentation in AMD subjects was proposed in [23]. In this method, a deformable model for accurate shape descriptions of fluid-filled regions was presented followed by a nonlinear anisotropic diffusion filter to decrease the effect of speckle noise. This method

was evaluated on 7 AMD subjects and then quantitative and qualitative (good, fair and poor extraction) analysis were carried out. Note that all following methods are considered as automated methods. In [3], a supervised method was proposed to identify fluidfilled regions and retinal layers on 110 Bscans of 10 patients with DME. In this method, fluid and retinal layer positions were estimated by a kernel regression-based classification approach. The classification method is then used for more accurate segmentation of retinal layer boundaries using a graph theory and dynamic programming framework. In [24], intraretinal and sub-retinal fluid regions were detected in 78 SDOCT volumes from 23 AMD patients by a supervised 3D method. This method computes the local differences from normal appearance by 23 extracted features in each layer locally and the normal ranges of layer-specific feature variations derived from 13 SD-OCT volumes of normal subjects. Fluid-associated abnormalities in OCT images of 15 OCT volumes from 15 AMD subjects undergoing intravitreal antiVEGF injection treatment were segmented by a supervised 3D method [25]. This method includes two main steps. Layer segmentation, candidate fluid identification and retinal OCT flattening were carried out in the first step. In the second step, a probability constrained graph search graph cut method was proposed to refine the candidate fluid regions. In [26], fluid filled boundaries on B-scans (x and y-axes) and C-scans (zaxis) were identified in DME subjects by an unsupervised method based on fuzzy level set. These boundaries were combined to generate 3D segmentation of retinal fluid. Then, morphological characteristics were used for the elimination of artefactual fluid regions. Identification of vascular shadowing was also used for removing such artefactuals. The accuracy of this method in retinal fluid segmentation was evaluated on 10 DME subjects. In [27], an automated method based on artificial neural network combined with a segmentation framework based on geodesic graph cut for retinal fluid segmentation from OCT images of AMD subjects was presented. This method was evaluated on 30 OCT volumes from 10 AMD subjects at 3 different treatment stages. In [28], a supervised 2D segmentation of cyst regions in AMD patients was presented based on k-means cluster analysis and k-nearest neighbor classifier. 31 volume scans collected during a 4-year period from one AMD patient with a serous retinal detachment were used to evaluate this method. Localization of cysts in OCT images of DME subjects was proposed using an unsupervised 2D method in [29]. In this approach, sub-retinal layers segmentation, dark regions detection, and solidity analysis steps were carried out for fluid segmentation. In [30], microcystic macular edema pixels were classified using a supervised 2D method based on random forest classifier with intensity and spatial features. Automatic segmentation of intraretinal cysts from 3D OCT images of DME patients was presented using an unsupervised method with three-dimensional curvelet transform [31]. This method was evaluated on 30 OCT volumes of 10 subjects at 3 different stages of treatment. In [32], an unsupervised 2D retinal cyst segmentation method was proposed based on bilateral filter for DME subjects. This approach was tested on 16 vitreoretinal patients and three control subjects. Detection of

cysts among DME subjects with non-linear complex diffusion and mathematical morphology operations was proposed by an unsupervised 2D method in [33]. Detection of pseudo-cysts on Spectralis OCT and Cirrus OCT scans of DME subjects was proposed based on a supervised 2D method in [34], [35]. These methods were evaluated on a small sample size of five subjects. Supervised 3D cystoid macular edema segmentation of macular hole images of DME patients was presented in [36]. This method was tested in 3D OCT images from 18 subjects with cystoid macular edema and macular hole. Finally, in [37], a higher-order constraint-based supervised 3D method was proposed based on label propagation for fluid-associated regions segmentation. Fluid/cyst regions can be revealed in AMD and DME subjects. Therefore, proposed methods for fluid/cyst segmentation are categorized as AMD and DME methods. Also, segmentation methods can be applied to 2D OCT slices and 3D OCT volumes which are referred as 2D and 3D methods, respectively. Furthermore, segmentation methods that use fluid/cyst regions labeled by experts are referred as supervised methods. Unsupervised methods segment fluid/cyst regions directly without requiring manually segmented regions. Finally, if user interaction is used in segmentation process the method is semi-automated; otherwise it is automated. A brief explanation of how our proposed method is similar to and differs from the related works is summarized as follows. The proposed method defines a new cost function in NS domain. Then, the cost function is minimized which leads to an unsupervised segmentation of 2D OCT slices of DME patients. Therefore, the proposed method can be considered in the same category of methods proposed in [29], [32], [33]. Similar to the proposed method, methods in [3], [30], [34], [35], [38] are 2D methods applicable to DME subjects although these methods use supervised procedures in segmentation process. Methods in [26], [31] address the same problem as the proposed method except these methods have been applied to 3D OCT volumes rather than 2D slices. [36], [37] also present DME methods which have been applied to 3D OCT volumes but in a supervised manner. Apart from the proposed methods for DME subjects, methods in [24], [25], [27], [39] have been proposed for AMD subjects using a supervised 3D procedure. Correlation between initial vision and vision improvement with automatically calculated retinal cyst volume in treated DME subjects was analysed in [40].

III. R EVIEW OF N EUTROSOPHIC I MAGES A. Neutrosophic Set Consider that X is a universal set in the neutrosophic domain and a set A is included in X. Each member x in A is described with three real subsets of [0, 1] named as T , I and F . Element x in set A is expressed as x(t, z, f ), where t, z and f vary in T , I and F , respectively. x(t, z, f ) could be interpreted as it is t% true, z% indeterminate, and f % false that x belongs to A. T , I and F could be considered as membership sets [7].

B. Neutrosophic Image For using the concept of NS in image processing, an image should be transformed into the neutrosophic domain. Although the general method for this transformation was proposed by Guo et al. [8], the method of transformation is completely dependent on the image processing application. In Section IV, we propose our transformation method that is appropriate for OCT segmentation. Consider an image g with L gray levels. g can be mapped into T , I and F sets. Thus, the pixel p(i, j) in g is transformed into P N S(i, j) = T (i, j), I(i, j), F (i, j)) or P N S(t, z, f ) in neutrosophic domain. T , I and F are considered as white, noise and black pixel sets, respectively. P N S(t, z, f ) means that this pixel is %t true to be a white pixel, %z to be a noisy pixel and %f to be a black pixel. T , I and F are computed as follows [8], [9]: g(i, j) − gmin (1) T (i, j) = gmax − gmin F (i, j) = 1 − T (i, j); (2) δ(i, j) − δmin I(i, j) = ; (3) δmax − δmin w

2 2 X 1 X g(i + m, j + n) g(i, j) = 2 w w w

(4)

m=− 2 n=− 2

δ(i, j) = |g(i, j) − g(i, j)| (5) where g is the gray scale image, g is a filtered version of the image g filtered with an averaging filter, w is the window size for the averaging filter, gmax and gmin are the maximum and minimum of the g, respectively, δ is the absolute difference between g and g, δmax and δmin are also the maximum and minimum values of δ, respectively. IV. P ROPOSED METHOD In this section, segmentation approach is proposed for fluid/cyst segmentation. A. Transform the images to neutrosophic domain and denoising For fluid/cyst segmentation, first the OCT image is transformed into NS by the proposed method. This method not only transforms the image into NS but also changes the gray level of the noisy pixels. In NS, indeterminacy is one of the important concepts, and in the image processing domain, it is interpreted as noise. Here, a new definition of indeterminacy is proposed. In the basic NS-based image segmentation method, indeterminacy set was defined by (3). In this definition, the greater the difference between each pixel with the mean of its neighbors in the square window, the greater the indeterminacy. In (3), higher indeterminacy is assigned to pixels near the OCT layers and the boundaries while these are not noisy pixels. For example, in Fig. 1, pixel P 1 is not a noisy pixel; however, a high indeterminacy value is assigned to this pixel because of its difference with the mean of its surrounding pixels. We propose a new definition of indeterminacy of the image g in Algorithm 1.

Algorithm 1 Indeterminacy computation 1: 2: 3: 4:

Consider a rectangular Gaussian filter with the dimension of [g1, g2]. Rotate the filter in 10 different directions to cover 180 degrees of rotation. Apply all filters to image T to compute 10 filtered images: F Ik=1...10 . Compute I as: I(i, j) = mink |T (i, j) âˆ’ F Ik (i, j)|

Fig. 2: Transformation to NS domain. (a): input OCT Bscan, (b): subset T, (c): subset F and (d): subset I. transformed OCT Bscan to NS domain. The obtained T and I sets in Algorithm 2 are used in the subsequent steps of segmentation. B. ROI segmentation with shortest path graph

Fig. 1: Illustration of the proposed indeterminacy definition. Instead of considering the difference between each pixel and the mean of its surrounding window, the minimum difference between each pixel and the mean of its neighbors in 10 different directions is considered. By this way, for pixels like P 2 (which has similar situation as P 1) in Fig. 1, the filter d1 (in horizontal direction) results in the lowest difference with P 2. Thus, the indeterminacy of this pixel is not increased. Pixels with high indeterminacy are considered as too noisy. A Îť-correction operation is proposed to decrease the noise effect as defined in (6)-(8). In this definition, the very noisy pixels are blurred with the filter which has the greatest difference (the biggest penalty is considered for these pixels). ( T (i, j), if I(i, j) < Îť TĚ„ (Îť) = (6) TÎť0 (i, j), otherwise Ind(i, j) = argmaxk |g(i, j) âˆ’ F Ik (i, j)|

(7)

TÎť0 (i, j) = F IInd(i,j) (i, j) (8) where F I has been defined in Algorithm 1. Based on this definition for indeterminacy set and Îťcorrection operation, the image g is transformed and denoised as described in Algorithm 2. Algorithm 2 Proposed image transformation to NS domain 1: 2:

3: 4: 5: 6:

Inputs: g (input image), Output: T , I and F . g(i, j) âˆ’ gmin Compute T (i, j) = 1 âˆ’ , F = 1 âˆ’ T and I gmax âˆ’ gmin with Algorithm 1. Apply Îť-correction operation to T set. Compute T = TĚ„ (Îť), F = 1 âˆ’ T and I with Algorithm 1. If |Entropy(Ik ) âˆ’ Entropy(Ikâˆ’1 )|< go to 6, otherwise go to 3. End.

The P entropy is defined as Entropy(I) = P m n âˆ’ i=1 j=1 I(i, j)log2 I(i, j). Since fluid regions are darker than other regions, the inverse of intensity is considered as T set. Therefore, pixels with high T (i, j) values are likely to be in fluid regions. Fig. 2 shows a

Fluid/cyst regions are all located between ILM and RPE layers in the OCT images of DME subjects. Therefore, the first step of our proposed method is the segmentation of these layers as the ROI in the NS domain. This step is very important due to two aspects. First, the background region is very similar to fluid/cyst regions in both brightness and texture. This can easily mislead the segmentation method since this method is based on an unsupervised clustering scheme and will be affected by the similarity between the desired fluid/cyst regions and the irrelevant background region. The second reason for ROI segmentation is speeding up since the ROI is processed instead of the whole image. The proposed ROI segmentation is derived from the method presented in Chiu et al. [19]; however, the method in [19] was proposed for the OCT layers of normal adult eyes while our proposed method is proposed for abnormal OCT images. For this task, the graph is constructed from each OCT image by mapping each pixel in the image to one node in a graph. We only consider the local relationship between pixels. Therefore, by considering the local relationship for 8 neighbors of each pixel, the 8-regular graph is constructed. For ILM segmentation, the image is first filtered with filter H for the calculation of vertical gradient of each pixel using (9). ďŁŽ ďŁš âˆ’2 V erGrad = T âˆ— H, H = ďŁ° 0 ďŁť (9) 2 where T is the transformed image in NS. The proposed weight computation between any two arbitrary pixels (a1 , b1 ) and (a2 , b2 ) is defined by (10): W ((a1 , b1 ), (a2 , b2 )) = 4 Ă— M axG âˆ’ V erGrad(a1 , b1 ) âˆ’ V erGrad(a2 , b2 ) + 2 Ă— mean(R) (10) where M axG is the maximum gray level of the image and R is a set of D pixels above (a1 , b1 ). In all experiments D is set to 40. Based on the filter H, the maximum of V erGrad is 2 Ă— M axG. So, the maximum of V erGrad(a1 , b1 ) + V erGrad(a2 , b2 ) is 4 Ă— M axG. Pixels which are located in the first layer have the maximum vertical gradient. Therefore, the minimum weight will be assigned between them and then they have the highest chance to be selected by the shortest path algorithm. The main problem is that there are pixels in other layers which also have the maximum vertical gradient.

Fig. 3: Weight computation parameters in (a): ILM, (b): ISM and (c): OPL segmentation.

where T and I represent true and indeterminacy sets in NS, respectively, which were computed in Algorithm 2, and Mijk and Nij represent the membership degree of pixel (i, j) to main cluster k and NC, respectively, K is the number of main clusters, Ck is the center of kth main cluster, m is a constant and n1 and n2 are the image dimension. The proposed cost function (13) is derived from fuzzy c-means clustering whose convergence has been proved in [41]. Jm =

N X K X

m Mij ||Xi − Cj ||2

(13)

i=1 j=1

Fig. 4: ROI segmentation for a sample Bscan. Fig. 3 (a) shows the pixels which have the maximum vertical gradient (and thus minimum weight) but not located in the first layer. The weight between these pixels is increased by introducing the term 2 ∗ mean(R); note that this term is very small for the pixels in ILM. Therefore, it is guaranteed that the pixels in ILM always have the minimum weight. After graph construction, the gradients of all boundary pixels are set to M axG and then the start and end points are selected randomly from nodes corresponding to the pixels in the first and the last columns of the image, respectively. The proposed weight computation equation assigns non-negative weights between nodes. Therefore, the Dijkstra shortest path algorithm can be applied to this problem to find the shortest path between start and end points and this path is considered as the ILM layer. The procedure for the segmentation of RPE is same except for the filter H which is inversed. Also, note that R corresponds to a set of D pixels under (a1, b1). A result of ROI segmentation for a sample OCT image is shown in Fig. 4. C. Segmentation scheme for fluid/cyst regions

where Xi and Cj represent data point i and cluster center j, respectively, and Mij represents the membership degree of data point i to cluster j. In fuzzy c-means, for each data point Xi , the highest membership is assigned to the closest cluster since the goal is to minimize the cost function. Our cost function includes two terms, one for the main clusters and another for NC. We consider two conditions for pixel (i, j) to have the highest membership degree to the main cluster k: (a) pixel (i, j) should have the minimum distance from the main cluster center k rather than other clusters, (b) pixel (i, j) should have the minimum indeterminacy. Similarly, there are also two conditions for pixel (i, j) to have the highest membership degree to NC:P (a) having the maximum sum distance from all K main clusters k=1 kTij − Ck k2 and (b) having the maximum indeterminacy. The maximum distance (in intensity) between any two pixels is 1 since all sets in NS have been normalized to PKthe interval [0,21]. Therefore, the maximum quantity for k=1 kTij − Ck k is K. For considering the constraint in (12), the Lagrange cost function in (14) is constructed. n1 X n2 X K X J(M, N, C) = (Iij Mijk )m kTij − Ck k2 i=1 j=1 k=1

In this research, an unsupervised clustering method is proposed for the segmentation of fluid/cyst regions in OCT images. For this task, a new cost function is proposed. Although the OCT images were denoised in Algorithm 2, in the clustering scheme, the effect of noise is still considered in the proposed cost function in two steps. In the first one, an extra cluster referred as noise cluster (NC) is considered beside K main clusters. Therefore, it is expected that the noisy pixels will be assigned to the NC. The second step is that the indeterminacy of the pixels is included in the cost function. The proposed cost function is presented in (11). J(M, N, C) n1 X n2 X K X = (Iij Mijk )m kTij

n1 X n2 X

K X

i=1 j=1

k=1

((1 − Iij )Nij )m ||K −

λij

||Tij − Ck ||2 ||−

(14)

n2 n1 X X (Mijk + Nij − 1) i=1 j=1

For minimizing this cost function, gradient descent approach is used. Therefore, 1 −2 m−1 λ ij m − 1 m ) ||Tij − Ck || m − 1 Mijk = Iij ( (15) m Nij = (1−Iij )

m−1 m

1 −2 K X λij m − 1 2 m−1 ( ) ||K − ||Tij −Ck || || m k=1

(16)

i=1 j=1 k=1

Pn1 Pn2

n1 X n2 K X X − Ck k2 + ((1 − Iij )Nij )m ||K − ||Tij − Ck ||2 || i=1 j=1

k=1

(11) s.t.

K X k =1

Mijk + Nij = 1

(12)

Ck =

m m i=1 j=1 (Iij Mijk ) ((1 − Iij )Nij ) Tij Pn1 Pn2 m m i=1 j=1 (Iij Mijk ) ((1 − Iij )Nij )

(17)

The final proposed clustering algorithm is described in Algorithm 3. The result of clustering Algorithm 3 for a sample Bscan by considering just fluid regions in step 5 is depicted in Fig. 5.

Algorithm 3 Proposed clustering algorithm 1: 2: 3: 4: 5:

Initialize K, m, Mijk and Nij . Compute Ck , Mijk and Nij using (15)-(17). If |Jk (M, N, C) âˆ’ Jkâˆ’1 (M, N, C)|< go to 4 else go to 2. Assign the cluster number with the highest membership to each pixel, the cluster number for NC is K + 1. Based on the definition of T in Algorithm 2, fluid pixels have the maximum value. Therefore, these pixels are in the cluster Ck with the greatest center value.

Fig. 5: Clustering result for a sample Bscan. D. Post processing with middle layers segmentation The proposed post-processing is an important step in image segmentation. After segmentation with the proposed clustering scheme, some segmented regions are incorrectly marked as fluid/cyst. The regions between OPL and ISM are similar to fluid/cyst regions in both texture and brightness. In the proposed approach, another layer segmentation method for OPL and ISM segmentation is proposed which is similar to the proposed ROI segmentation method and is based on the Dijkstra graph shortest path algorithm. After the segmentation of fluid regions, the regions between OPL and ISM layers are ignored during post-processing. For the segmentation of ISM, the graph is constructed from the image similar to that for ROI, except for the weight computation between nodes which is proposed to be calculated by (18). W (p1 , p2 ) = 2 Ă— M axG âˆ’ [mean(L1) âˆ’ mean(U 1)] (18) âˆ’ [mean(L2) âˆ’ mean(U 2)] where L1, U 1, L2 and U 2 are the set of pixels under p1 , above p1 , under p2 and above p2 , respectively, and are shown in Fig. 3(b). For OPL segmentation, the weights are computed using (19). W (p1, p2) = 2 Ă— M axG âˆ’ [mean(U 1) âˆ’ mean(L1)] (19) âˆ’ [mean(U 2) âˆ’ mean(L2)] + Î˛ Ă— D. where D is the distance of p1 from ISM layer, U1 , L1 , U2 and L2 are same with (18). Fig. 3(c) shows all parameters for weight computation in OPL. The extra term Î˛ Ă— D in (19) increases the weight of pixels that are far from ISM. Therefore, this term enforces the shortest path algorithm to be as close as possible to ISM. This step is necessary since all the regions between ISM and OPL are ignored in the post processing step. Therefore, the proposed method tries to find OPL layer under the fluid regions. For both OPL and ISM, the start and end points for finding shortest path are selected from any arbitrary pixels in the first and the last columns of ROI image. Finally, OPL is flattened using (20).

Fig. 6: OPL and ISM segmentation result for a sample Bscan. ( Oi , ifabs(Oi âˆ’ Ii ) <= mean(Dif ) Oi = (20) Ii + AD, otherwise where Oi and Ii represent OPL and ISM pixels in ith column, respectively, Dif is a vector that shows the distance between OPL and ISM in each column and AD is the average distance between OPL and ISM. OPL and ISM segmentation results are shown in Fig. 6. E. Fluid/cyst segmentation 1) Automated determination of the number of clusters (K): In the proposed clustering method, the number of main clusters (K) affects the segmentation results significantly. Therefore, for having appropriate segmentation results, K should be determined carefully. The main question is what K is appropriate for any arbitrary OCT scan. For this task, all ground truth images were analyzed and we found that for images with a great deal of fluid/cyst, ophthalmologists are interested in the segmentation of the image with more general and bigger segments. Therefore, segmentation with a small K has the best correlation with the segmentation annotated by an ophthalmologist. For the images with smaller fluid/cyst regions, a larger K is more appropriate. By this interpretation, an automated method is proposed for determining K. In this method at first the minimum and maximum values for K are considered and then the fluid/cyst regions are estimated. Then, an appropriate K is determined based on this estimation. Algorithm 4 describes the proposed method for automated determination of K: Algorithm 4 Automated determination of K 1: 2: 3: 4: 5:

Find the mean (Âľ) and standard deviation (Ďƒ) of the normal distribution of the pixel intensities between OPL and ISM. Find the histogram of pixels between ILM and OPL regions. Count the number of pixels in the intensity interval [Âľ âˆ’ Ďƒ, Âľ + Ďƒ] and name it as A. Denote the maximum and minimum of A for all OCT images in each subject as Amin and Amax . Map the A of the each OCT image to appropriate K A âˆ’ Amin ,K = with the following Eqs: A1 = 1 âˆ’ Amax âˆ’ Amin round((Kmax âˆ’ Kmin )A1 + Kmin )

2) Final fluid/cyst segmentation algorithm: Based on proposed methods in previous sections, the final fluid/cyst segmentation is performed using the following steps: 1) Transform OCT scans to NS domain using Algorithms 1 and 2. 2) Segment ILM and RPE as ROI and ISM and OPL as middle layers. 3) Compute appropriate number of clusters using

Fig. 7: Final segmentation result for a sample Bscan. Algorithm 4. 4) Segment fluid/cyst regions using Algorithm 3. 5) Remove the regions between ISM and OPL and small regions under H pixles. Final segmentation result for a sample Bscan is shown in Fig. 7. V. E XPERIMENTAL SETUP A. Datasets In this work, three datasets have been used for the evaluation of the proposed algorithms. The first dataset is a local dataset from the UMN ophthalmology clinic containing 725 images from 29 DME subjects which were taken using the Heidelberg Spectralis imaging system. Each image is obtained by averaging 12-19 frames with the resolution of 5.88Âľm/pixel along the length and 3.87Âľm/pixel along the width. Fluid/cyst regions were segmented by two UMN ophthalmologists (DDK and PMD). The second dataset from Duke includes ten DME subjects (patients) with lateral and azimuthal resolutions ranging from 10.94 to 11.98Âľm/pixel and 118 to 128Âľm/pixel, respectively [3]. This dataset is available online and includes automated and manual segmentation results1 . The third dataset is from the OPTIMA Cyst Segmentation Challenge and contains 4 DME subjects with 49 images per subject where the image resolution varies from 512x496 to 512x10242 . B. Parameters tuning Cluster numbers (K) are computed adaptively by Algorithm 4. In Algorithm 1: parameters [3, 9] were used for filter dimension [g1, g2]. These dimensions were selected to just have a rectangular filter. = 0.001 in Algorithm 2 and = 0.01 in Algorithm 3 were selected in such small quantities to make sure that there is no significant change in the cost function J and entropy of the indeterminacy set I. Parameters Îť = 0.7 in Algorithm 2 and Î˛ = 0.008 in layer segmentation were determined by experiments. Parameters m, M and N in Algorithm 3 are initialized randomly with the same procedure as in conventional fuzzy c-means algorithms. In Algorithm 4, the minimum and maximum number of clusters were set to 4 and 10, respectively. The reason is that in our experiments, if these numbers are selected to be smaller than 4 and greater than 10, general (big segments) and very small fluid/cyst regions are segmented, respectively, which do not correlate well with manual expert segmentations. In layer segmentation methods, R was set to 40. This quantity includes while pixels in the rectangular area above the pixels in the ISM layer. Parameters U1 , L1 , U2 and L2 were set to 10 to only consider dark pixels above ISM and bright pixels below ISM. The difference between the means of these two pixel sets is high and it reduces the weights of the pixels located in ISM 1 http://duke.edu/

sf59/ChiuË™BOEË™2014Ë™dataset.htm

2 http://optima.meduniwien.ac.at/challenges/optima-segmentation-

challenge-1/

Fig. 8: Final fluid/cyst segmentation results in Duke dataset. (a) input OCT images, (b): expert 1 segmentation, (c): expert 2 segmentation, (d): method in [3] and (e): proposed method segmentation.

Fig. 9: Final fluid/cyst segmentation results in Optima dataset. (a) input OCT images, (b): expert 1 segmentation, (c): expert 2 segmentation and (d): proposed method segmentation. layer by equation (20). It means that pixels in ISM layer are good candidates to be selected by the shortest path algorithm. However the layer segmentation method is not very sensitive to the small changes of these parameters since the means of these pixel sets are used in weight computation equations. For example, if we assign 45 to R instead of 40, the mean of 45 pixels is not very different from the mean of 40 pixels. Finally, in post-processing step, H was set to 50. It was selected due to the fact that in our datasets, ophthalmologists did not segment very small regions under 50 pixels in area as fluid/cyst region. Therefore, all segmented regions under 50 pixels are ignored in post-processing. VI. E XPERIMENTAL RESULTS Our proposed fluid/cyst segmentation has been tested on the three mentioned datasets. The final fluid/cyst segmentation results on Duke, Optima and local UMN datasets are illustrated in Figs. 8, 9 and 10, respectively. The proposed fluid/cyst segmentation method has been evaluated with respect to dice coefficient, precision and sensitivity criteria. Accuracy and specificity criteria havenâ€™t been used

Fig. 10: Final fluid/cyst segmentation results in UMN dataset. (a) input OCT images, (b): expert 1 segmentation, (c): expert 2 segmentation and (d): proposed method segmentation. here since these are biased to very high values (close to 100%). This is due to the large number of negative (non-fluid) pixels and the large number of true negative (TN) segmented pixels. Segmentation results of our proposed method using Duke dataset are compared with manual segmentation results by two ophthalmologists and with the results of the method in [3]. The dice coefficients, sensitivity and precision for all subjects in each dataset are reported in Tables I-III. Our proposed method achieves an average dice coefficient of 57.51% and outperforms the method in [3] that achieves an average dice coefficient of 49.24%. Furthermore, our proposed method achieves the highest dice coefficients of 70.52% for the Optima dataset. The method in [3] achieves only about 1% better sensitivity than our method while the proposed method achieves about 8% and 5% higher dice coefficient and precision, respectively, compared to [3]. It should be noted that the method in [3] achieved an average dice coefficient of 78% for fluid, NFL, IPL, INL, OPL, ISM, ISE and RPE segmentation. However, we compared our segmentation results with only fluid segmentation results of [3] since layer segmentation is not the main contribution of this paper and is just used for pre-processing and post-processing steps. The goal of the layer segmentation method is to avoid false positives in fluid/cyst segmentation, and not to segment layers that correlate well with annotations by the ophthalmologist. Therefore, segmented sub-retinal layers have different positions as compared with segmented layers by the ophthalmologist and these are not comparable. Note that the approach in [3] for Duke dataset required 11.4s compared to 29.25s for the proposed approach. However, this segmentation is significantly less than 5.5 minutes required for manual segmentation by an expert. Also, our segmentation results in Optima dataset are compared with the method in [31] (which achieved the third rank in Optima challenge) in TABLE III. Note that the first and second ranks have not been published yet. Our proposed method achieves 70.52%, 88.84% and 73.89% with respect to dice coefficient, sensitivity and precision, respectively, which outperforms the method in [31] with 63.63%, 66.49% and 50.03%. Finally, the proposed method achieves 69.40%, 76.79% and 74.91% with respect to dice coefficient, sensitivity and precision, respectively, in the UMN dataset. The results of different steps in the proposed fluid/cyst segmentation scheme for sample Bscans and the average run times of the different steps of the segmentation algorithm have been reported in supplementary parts of the

Fig. 11: Averge Dice Coefficients, Sensitivity and Precision of methods in all datasets. TABLE I: Dice Coefficients, sensitivity and precision of all subjects in local collected UMN dataset. Dice Coefficients Exp. 1 Exp. 2 83.32% 79.25% 83.58% 83.50% 57.24% 49.81% 37.01% 37.07% 67.45% 67.45% 83.45% 82.88% 75.80% 75.80% 66.76% 62.74% 65.65% 65.21% 92.00% 88.00% 78.23% 73.72% 60.81% 61.98% 62.39% 62.42% 65.05% 65.00% 67.31% 67.77% 66.65% 66.24% 66.29% 66.29% 64.31% 64.66% 81.37% 81.29% 66.24% 66.16% 73.20% 73.08% 67.53% 67.52% 76.83% 76.64% 76.24% 76.24% 73.65% 72.95% 70.00% 70.03% 59.96% 58.04% 69.37% 69.77% 67.89% 68.02% 69.85% 68.95%

S1 S2 S3 S4 S5 S6 S7 S8 S9 S 10 S 11 S 12 S 13 S 14 S 15 S 16 S 17 S 18 S 19 S 20 S 21 S 22 S 23 S 24 S 25 S 26 S 27 S 28 S 29 Ave.

Sensitivity Exp. 1 Exp. 2 94.08% 89.68% 86.76% 86.06% 85.84% 74.40% 50.26% 49.61% 78.62% 78.62% 85.06% 83.88% 86.99% 86.99% 79.44% 75.38% 76.55% 75.97% 100.00% 96.00% 79.49% 74.14% 57.96% 58.73% 82.68% 82.72% 82.36% 82.06% 70.77% 67.89% 69.78% 69.19% 67.46% 67.46% 72.77% 72.57% 83.33% 82.95% 68.57% 68.31% 66.88% 66.75% 66.29% 66.28% 84.38% 83.92% 82.87% 82.87% 74.01% 72.25% 82.36% 82.55% 58.62% 58.77% 88.13% 86.97% 85.03% 83.86% 77.49% 76.10%

Precision Exp. 1 Exp. 2 83.14% 83.21% 85.39% 85.52% 62.34% 62.32% 50.10% 50.41% 65.68% 65.68% 84.16% 84.25% 78.63% 78.63% 73.56% 73.64% 80.65% 80.52% 92.00% 92.00% 89.69% 90.07% 82.39% 82.91% 61.60% 61.59% 62.31% 62.41% 71.03% 72.25% 65.67% 65.50% 73.69% 73.69% 66.28% 67.13% 82.48% 82.52% 76.74% 76.92% 93.17% 93.17% 78.54% 78.54% 80.15% 80.27% 76.04% 76.04% 75.80% 76.44% 71.37% 71.03% 65.21% 61.23% 73.58% 76.14% 69.04% 70.34% 74.84% 74.98%

paper1 . VII. S TATISTICAL C OMPARISON In the experimental results section, the average of each measure was used as an overall performance criterion. In this section a statistical test named as pairwise U-MannWhitney test is used for evaluating which method has better performance statistically [42]. Different pairwise comparisons between our proposed method and [3] using different criteria are carried out using the mentioned test. The confidence level threshold of 90% is considered. Test results show that the 1 The

Address of the supplementary materials.

TABLE II: Dice Coefficients, sensitivity and precision of all subjects in Duke dataset.

S1 S2 S3 S4 S5 S6 S7 S8 S9 S 10 Ave.

Dice Cofficients Method in [3] Prop. Method Exp. 1 Exp. 2 Exp. 1 Exp. 2 61.79% 42.38% 74.53% 52.96% 48.10% 43.10% 62.56% 59.38% 52.98% 37.05% 72.90% 42.69% 42.81% 42.69% 42.87% 43.90% 66.19% 57.24% 61.49% 60.64% 49.67% 57.43% 48.81% 55.86% 52.38% 45.84% 80.56% 74.01% 48.87% 57.65% 55.25% 59.35% 51.09% 14.78% 58.33% 20.78% 59.40% 53.34% 66.27% 56.92% 53.33% 45.15% 62.36% 52.65%

Sensitivity Method in [3] Prop. Method Exp. 1 Exp. 2 Exp. 1 Exp. 2 68.03% 48.47% 79.28% 55.40% 74.86% 70.33% 70.91% 68.76% 92.99% 73.15% 83.63% 53.96% 79.76% 78.32% 70.86% 74.05% 81.09% 71.47% 81.58% 72.62% 79.03% 77.92% 64.85% 61.73% 85.79% 87.76% 89.72% 94.52% 60.48% 61.41% 59.07% 55.34% 56.07% 16.36% 58.29% 18.40% 54.47% 45.69% 73.75% 58.64% 73.26% 63.09% 73.19% 61.34%

Precision Method in [3] Prop. Method Exp. 1 Exp. 2 Exp. 1 Exp. 2 57.91% 58.66% 71.77% 70.99% 44.30% 40.83% 65.79% 60.21% 47.49% 50.08% 72.18% 71.69% 34.93% 38.21% 39.33% 43.61% 60.42% 62.74% 56.05% 58.74% 48.25% 59.53% 55.93% 69.64% 58.96% 54.46% 82.39% 74.12% 64.74% 75.87% 65.73% 76.14% 84.87% 88.36% 66.91% 71.05% 72.68% 84.67% 61.01% 66.32% 57.46% 61.34% 63.71% 66.25%

TABLE III: Dice Coefficients, sensitivity and precision of all subjects in Optima dataset.

S1 S2 S3 S4 Ave.

Dice Coefficients Exp. 1 Exp. 2 83.38% 81.86% 59.53% 57.46% 71.34% 79.09% 65.90% 65.59% 70.04% 71.00%

Sensitivity Exp. 1 87.59% 97.17% 85.39% 87.19% 89.34%

Exp. 2 83.03% 92.56% 91.57% 86.24% 88.35%

Precision Exp. 1 84.45% 59.51% 77.07% 71.62% 73.16%

Exp. 2 87.30% 59.61% 79.21% 72.36% 74.62%

proposed method with the mean rank of 12.9 statistically outperforms the method in [3] with the mean rank of 8.9 in terms of dice coefficient with a confidence level of 93%. For sensitivity and precision all confidence levels are under threshold; therefore, which method is significantly better cannot be concluded, although, for example, in precision measure the mean rank of the proposed method and the method in [3] are 12.1 and 8.9, respectively, with a confidence level of 77.4%. VIII.

DISCUSSION

A. The effect of hard exudates and signal intensity In some DME eyes, a collection of hard exudates exists which are imaged as intensely hyper-reflective structures. With regard to the evaluation of the proposed methods in subjects with these structures, we collected 2 subjects with these structures from the UMN eye clinic. In Fig. 12, (a) shows a sample of an OCT image with hyper-reflective structures and (b) is the transformed image corresponding to subset T in neutrosophic domain. Based on the step 2 of Algorithm 2 and the result of transformation in Fig. 12(b), very low memberships are assigned to hyper-reflective regions. Also, in indeterministic set, higher memberships are assigned to these regions (see Fig. 12(c)) which means that in the designed cost function in the clustering scheme, these are considered as noisy pixels and assigned to noise cluster. Therefore, these structures do not affect the fluid/cyst segmentation results. With regard to the effect of hard exudates on the proposed layer segmentation algorithms, although these regions create high vertical-gradient points (good candidates for being selected by shortest path algorithms), the shortest path does not include these regions. The cost of a path through these regions would be much higher as it would contain many low vertical-gradient nodes. Fig. 12(d) shows the segmented OPL layer and fluid/cyst regions by proposed algorithms.

Fig. 12: Hyper-reflective structures in a sample OCT Bscan. (a): Input OCT, (b): transformed image to neutrosophic domain (c): indeterminacy set and (d): segmented OPL and fluid/cyst regions.

Fig. 13: Signal intensity effect. It is clear that the intensity of OCT Bscans may be different in even subsequent slices of an OCT volume. In our collected dataset and other two datasets, there are variations in signal intensities. In the proposed Algorithm 2 for transforming Bscans to neutrosophic domain, the highest membership is assigned to fluid/cyst regions. This membership is computed based on the intensity of each individual Bscan. In the proposed clustering Algorithm 3, the cluster with the highest center value in T subset in neutrosophic domain is designed for fluid/cyst regions. It guarantees that if the Bscan be either darker or brighter in subsequent slices, fluid regions always are located in the cluster with maximum membership. Fig. 13 shows the two subsequent slices in one of the subjects of Duke dataset. It can be observed that the Bscan in right side with lower noise level in background is darker than left side. Segmented fluid/cyst results for these samples show that the proposed method is robust with respect to the variation of signal intensity. In these examples, although the average of the points in fluid/cyst clusters in the T set are 0.71 and 0.84 for the left and right Bscans, respectively, it does not affect the fluid/cyst segmentation results. In contrast with semiautomated and supervised methods reported in section II which are sensitive to initial seed points and training samples for fluid regions (due to the variation of intensity level for fluid regions in different Bscans), our proposed method is fully-automated, unsupervised, more robust and does not need any user seed point or training set for fluid regions. It should be noted that in some cases with abnormal intensity changes in layer structures, our proposed layer segmentation method can lead to a part of the layer path that is misleading. Figs. 14 (a) and (b) show two examples of such

TABLE IV: The effect of NS in segmentation results.

UMN dataset

Duke dataset

Optima dataset

Fig. 14: Middle layer segmentation errors. cases. This is due to the existence of the high vertical-gradient pixels which stem from the intensity changes from very bright regions to darker ones. However, these errors do not affect the final fluid/cyst segmentation results since segmented layers are just used in the ROI segmentation and post processing steps. In few images which are very dim, gradient information of the layer structures is almost lost and the OPL layer segmentation method follows a wrong path. This path may be even under the ISM layer. Fig. 14 (c) shows an example of such case. In these cases, the final fluid/cyst segmentation results will be affected and false positive pixels are segmented as fluid/cyst region which is due to the fact that in very dim images many pixels are similar to fluid/regions and the middle layer segmentation fails to ignore them in the post-processing step.

Dice Coff. Sensitivity Precision Dice Coff. Sensitivity Precision Dice Coff. Sensitivity Precision

With NS 69.40% 76.79% 74.91% 57.51% 67.27% 64.98% 70.52% 88.84% 73.89%

Without NS 64.76% 70.45% 68.14% 48.34% 59.73% 58.15% 59.12% 73.17% 61.08%

Fig. 15: Transformation to NS domain. (a): input OCT Bscan, transformed image by: (b): proposed method, (c): method in [8] and (d): method in [9]. 11.14% improvement in dice coefficient which is due to the existence of higher amount of speckle noise in this dataset. The minimum improvement of segmentation accuracy is 4.64% in the UMN dataset which has the lowest level of noise. It can be concluded that using proposed transformation to NS domain has more benefits for images with lower signal to noise ratio.

B. Benefits of NS theory in OCT analysis Here, the effect of the proposed indeterminacy definition and λ-correction operation in NS domain is analyzed qualitatively and quantitatively. Although transforming image to NS domain is not a direct denoising approach, the effect of noise is decreased in this transformation by operations such as α-mean and β-enhancement in the conventional NS-based segmentation approaches and indeterminacy definition and λcorrection in the proposed method. For the evaluation of the proposed transformation method, the result of a transformed image by our method is compared with other transformation methods in [8], [9] which are based on NS. Fig. 15 (a), (b), (c) and (d) show the input OCT Bscan, transformed image by: our method, method in [8] and method in [9], respectively. It is visually clear that the sub-retinal layers structure is not preserved in [8]. Method in [9] preserves layer structure better than method in [8] although some artifacts are created near the hard exudates regions. Our proposed method preserves the layer structure without creating any artifact in the NS domain. It should be noted that we compare subset F in our method with the subsets T in [8], [9] since in Algorithm 2 subset T is inverted to assign highest memberships to pixels in fluid/cyst regions. Beside the benefits of using NS for hard exudates, hyperreflective structures and signal intensity variations, the effect of NS is evaluated quantitatively to show that how NS affects the final fluid/cyst segmentation results. TABLE IV shows how segmentation results are affected when Algorithms 3 and 4 are used in image domain for segmentation. It should be noted that in clustering Algorithm 3, subsets T is replaced with input image g and I is ignored. As it is reported in TABLE IV the best effect of using NS is observed in the Optima dataset with

C. Inter-observer variability analysis In this study, the automatic segmentation results are compared with two experts referred as grader 1 (G1) and grader 2 (G2). Correlation between observers is analyzed and reported in TABLE V. When the observer correlation is high, the segmentation errors are due to the segmentation algorithm, not inter-observer variation. The observer correlation of the used datasets in this study are 97.70%, 91.15% and 58.08% with respect to the dice coefficients in the UMN, Optima and Duke datasets, respectively. It can be concluded from TABLE V that the dice coefficient errors of our automated algorithm in UMN and Optima datasets are from the algorithm, not inter-observer variability while this error is increased by interobserver variability in the Duke dataset. D. The effect of the number of clusters in segmentation results The number of main clusters (K) affects the segmentation results significantly. The smaller K, the bigger clusters are created and vice versa for larger K. The segmentation results TABLE V: Inter-observer variability analysis.

UMN dataset

Duke dataset

Optima dataset

Dice Coff. Sensitivity Precision Dice Coff. Sensitivity Precision Dice Coff. Sensitivity Precision

Average of Automatic vs. G1 and G2 69.40% 76.79% 74.91% 57.51% 67.27% 64.98% 70.52% 88.84% 73.89%

Inter-Observer 97.70% 97.46% 99.34% 58.08% 60.20% 76.80% 91.15% 91.39% 94.35%

Automatic vs. G1 and G2 Intersection 69.88% 77.64% 74.62% 57.99% 81.68% 52.50% 71.80% 93.36% 70.93%

TABLE VII: Effect of removing the regions between OPL and ISM and small regions.

UMN Dataset

Fig. 16: Fluid/cyst segmentation results for (a) K=4 and (b) K=7.

Duke Dataset

Optima Dataset

TABLE VI: The effect of the number of clusters in segmentation results. Random Bscans Bscans with hard exudates Bscans with intensity changes Bscans without fluid

K=4 45.66% 33.73% 69.81% 21.13%

K=5 49.83% 31.59% 75.87% 25.19%

K=6 44.33% 44.68% 66.55% 26.49%

K=7 23.80% 49.52% 66.87% 48.84%

K=8 39.42% 38.76% 59.35% 59.42%

K=9 31.75% 35.33% 45.41% 75.51%

K=10 27.19% 41.46% 40.83% 88.64%

Automated K 64.69% 68.76% 86.64% 84.33%

for two different K values are depicted in Fig. 16. In this section the effect of K in final segmentation results is analyzed quantitatively. If K is considered as a constant value for all OCT scans, average segmentation accuracy is decreased. For example, for the scan shown in Fig. 16, K = 4 is more appropriate than K = 7. In scans with less fluid regions K = 7 is better. It means that for each scan, the corresponding K should be considered. To show that how the proposed method for the determination of K affects segmentation results, four groups of OCT scans (100 scans per group) including random scans, scans with hard exudates, scans with intensity changes and scans without fluid are considered. The main reason for this consideration is to evaluate the effect of K in all types of OCT scans (first group) and OCT scans which have more challenges for fluid segmentation (last three groups). The average dice coefficient of each group for different values of K (from the minimum to maximum values) is reported in TABLE VI. From TABLE VI it is concluded that for the first three groups, using the K obtained from the automated method improves segmentation results significantly. For the last group, considering the maximum value for K is better than automated K. This stems from the fact that segmentation with the maximum K segments very small fluid regions which are ignored in the last segmentation step in Section IV.E.2. E. Remove the regions between OPL and ISM and small regions To improve the fluid/cyst segmentation results, two heuristics are presented which ignore the incorrectly segmented regions as fluid/cyst. In the first one (heuristic 1), the regions between OPL and ISM are ignored. This idea stems from the fact that these regions are created from OPL elevation which is because of fluid/regions in the central part of OCT scans. These regions are similar to fluid/cyst regions in both texture and brightness and mislead the proposed segmentation method. In the second one (heuristic 2), small regions under H pixels are ignored. The reason for this heuristic is that the proposed segmentation method does not consider spatial dimension and connectivity of segmented regions. Therefore, small isolated regions are segmented as fluid/cyst which are false positive pixels. To show how these heuristics improve the segmentation results, TABLE VII reports the dice coefficient, sensitivity and precision of the proposed segmentation method

Dice Coff. Sensitivity Precision Dice Coff. Sensitivity Precision Dice Coff. Sensitivity Precision

Without post-processing

Ignore small regions

Ignore regions between OPL and ISM

61.79% 77.61% 71.32% 46.11% 67.95% 58.45% 62.19% 90.13% 70.56%

63.76% 77.13% 72.43% 50.65% 66.34% 61.78% 66.98% 90.66% 71.13%

68.33% 75.41% 72.99% 55.86% 67.38% 61.75% 65.72% 89.54% 70.33%

Ignore small regions and the regions between OPL and ISM 69.40% 76.79% 74.91% 57.51% 67.27% 64.98% 70.52% 88.84% 73.89%

with and without these heuristics separately. Both heuristics decrease the number of false positive pixels since the proposed segmentation method consider some regions between OPL and ISM and some small regions far from fluid/cyst regions as fluid/cyst (see Fig. 5 in this document and Fig. 1 (b) in supplementary materials). Therefore, as shown in TABLE VII, applying these heuristics improves dice coefficient by 9%, 11% and 8% in UMN, Duke and Optima datasets, respectively. Also, precisions are increased 3%, 6% and 3% in UMN, Duke and Optima datasets, respectively. With respect to sensitivity measure, there is no considerable improvement in any metric since the number of false negatives is not affected by such heuristics significantly. Finally, in both heuristics, in some cases the number of true positive pixels is slightly decreased which leads to a decrease in the sensitivity. IX. C ONCLUSION In this research, a fully automated algorithm for fluid/cyst segmentation in OCT images of the retina with DME pathology has been proposed based on graph shortest paths and neutrosophic sets. Accurate segmentation of DME biomarkers is important since it can provide a quantitative measure for diagnosis of DME. To show the efficiency of the proposed method, it was tested on three OCT datasets with DME. Segmentation results show that the segmented images obtained by our proposed algorithm are not only in close agreement with the manual segmentations of the two ophthalmologist experts but also achieve better performance with respect to dice coefficient and precision comparison criteria, as compared to a prior method [3]. Future efforts will be directed towards the segmentation of fluid-filled regions in age-related macular degeneration (AMD). Future efforts will also be directed towards fine-tuning the algorithm for OCT images obtained from other manufacturers. Finally, reproducibility studies between segmentation following repeat imaging can be addressed as another future work. X. DATA ACCESS The local UMN DME dataset and the segmented images of the approach presented in this paper for the UMN, Duke and Optima data sets are available from Prof. Keshab Parhiâ€™s web site: http://people.ece.umn.edu/âˆźparhi/DATA/OCT/DME/. ACKNOWLEDGMENT This research was supported in part by the Minnesota Lions and an unrestricted grant from Research to Prevent Blindness.

R EFERENCES [1] D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito et al., “Optical coherence tomography,” Science (New York, NY), vol. 254, no. 5035, p. 1178, 1991. [2] R. Kafieh, H. Rabbani, M. D. Abramoff, and M. Sonka, “Intra-retinal layer segmentation of 3d optical coherence tomography using coarse grained diffusion map,” Medical image analysis, vol. 17, no. 8, pp. 907– 928, 2013. [3] S. J. Chiu, M. J. Allingham, P. S. Mettu, S. W. Cousins, J. A. Izatt, and S. Farsiu, “Kernel regression based segmentation of optical coherence tomography images with diabetic macular edema,” Biomedical optics express, vol. 6, no. 4, pp. 1172–1194, 2015. [4] J. W. Yau, S. L. Rogers, R. Kawasaki, E. L. Lamoureux, J. W. Kowalski, T. Bek, S.-J. Chen, J. M. Dekker, A. Fletcher, and J. Grauslund, “Global prevalence and major risk factors of diabetic retinopathy,” Diabetes care, vol. 35, no. 3, pp. 556–564, 2012. [5] J. Pe’er, R. Folberg, A. Itin, H. Gnessin, I. Hemo, and E. Keshet, “Upregulated expression of vascular endothelial growth factor in proliferative diabetic retinopathy.” British journal of ophthalmology, vol. 80, no. 3, pp. 241–245, 1996. [6] D. C. DeBuc, A review of algorithms for segmentation of retinal image data using optical coherence tomography. INTECH Open Access Publisher, 2011. [7] F. Smarandache, A Unifying Field in Logics Neutrosophic Logic. Neutrosophy, Neutrosophic Set, Neutrosophic Probability. American Research Press, 2005. [8] Y. Guo and H.-D. Cheng, “New neutrosophic approach to image segmentation,” Pattern Recognition, vol. 42, no. 5, pp. 587–595, 2009. [9] M. Zhang, L. Zhang, and H. Cheng, “A neutrosophic approach to image segmentation based on watershed method,” Signal Processing, vol. 90, no. 5, pp. 1510–1517, 2010. [10] A. Sengur and Y. Guo, “Color texture image segmentation based on neutrosophic set and wavelet transformation,” Computer Vision and Image Understanding, vol. 115, no. 8, pp. 1134–1144, 2011. [11] A. Heshmati, M. Gholami, and A. Rashno, “Scheme for unsupervised colour–texture image segmentation using neutrosophic set and nonsubsampled contourlet transform,” IET Image Processing, vol. 10, no. 6, pp. 464–473, 2016. [12] Y. Guo, A. Şengür, and J. Ye, “A novel image thresholding algorithm based on neutrosophic similarity score,” Measurement, vol. 58, pp. 175– 186, 2014. [13] J. Shan, H. Cheng, and Y. Wang, “A novel segmentation method for breast ultrasound images based on neutrosophic l-means clustering,” Medical physics, vol. 39, no. 9, pp. 5669–5682, 2012. [14] Y. Guo and A. Şengür, “A novel image edge detection algorithm based on neutrosophic set,” Computers & Electrical Engineering, vol. 40, no. 8, pp. 3–25, 2014. [15] Y. Guo and A. Sengur, “Ncm: Neutrosophic c-means clustering algorithm,” Pattern Recognition, vol. 48, no. 8, pp. 2710–2724, 2015. [16] B. Peng, L. Zhang, and D. Zhang, “A survey of graph theoretical approaches to image segmentation,” Pattern Recognition, vol. 46, no. 3, pp. 1020–1038, 2013. [17] M. K. Garvin, M. D. Abràmoff, R. Kardon, S. R. Russell, X. Wu, and M. Sonka, “Intraretinal layer segmentation of macular optical coherence tomography images using optimal 3-d graph search,” IEEE transactions on medical imaging, vol. 27, no. 10, pp. 1495–1505, 2008. [18] Q. Yang, C. A. Reisman, Z. Wang, Y. Fukuma, M. Hangai, N. Yoshimura, A. Tomidokoro, M. Araie, A. S. Raza, D. C. Hood et al., “Automated layer segmentation of macular oct images using dual-scale gradient information,” Optics express, vol. 18, no. 20, pp. 21 293–21 307, 2010. [19] S. J. Chiu, X. T. Li, P. Nicholas, C. A. Toth, J. A. Izatt, and S. Farsiu, “Automatic segmentation of seven retinal layers in sdoct images congruent with expert manual segmentation,” Optics express, vol. 18, no. 18, pp. 19 413–19 428, 2010. [20] K. Lee, M. Niemeijer, M. K. Garvin, Y. H. Kwon, M. Sonka, and M. D. Abramoff, “Segmentation of the optic disc in 3-d oct scans of the optic nerve head,” IEEE transactions on medical imaging, vol. 29, no. 1, pp. 159–168, 2010. [21] S. J. Chiu, J. A. Izatt, R. V. O’Connell, K. P. Winter, C. A. Toth, and S. Farsiu, “Validated automatic segmentation of amd pathology including drusen and geographic atrophy in sd-oct images,” Investigative ophthalmology & visual science, vol. 53, no. 1, pp. 53–61, 2012.

[22] K. K. Parhi, A. Rashno, B. Nazari, S. Sadri, H. Rabbani, P. Drayna, and D. D. Koozekanani, “Automated fluid/cyst segmentation: A quantitative assessment of diabetic macular edema,” Investigative Ophthalmology & Visual Science, vol. 58, no. 8, pp. 4633–4633, 2017. [23] D. C. Fernandez, “Delineating fluid-filled region boundaries in optical coherence tomography images of the retina,” IEEE transactions on medical imaging, vol. 24, no. 8, pp. 929–945, 2005. [24] G. Quellec, K. Lee, M. Dolejsi, M. K. Garvin, M. D. Abramoff, and M. Sonka, “Three-dimensional analysis of retinal layer texture: identification of fluid-filled regions in sd-oct of the macula,” IEEE transactions on medical imaging, vol. 29, no. 6, pp. 1321–1330, 2010. [25] X. Chen, M. Niemeijer, L. Zhang, K. Lee, M. D. Abràmoff, and M. Sonka, “Three-dimensional segmentation of fluid-associated abnormalities in retinal oct: probability constrained graph-search-graph-cut,” IEEE transactions on medical imaging, vol. 31, no. 8, pp. 1521–1531, 2012. [26] J. Wang, M. Zhang, A. D. Pechauer, L. Liu, T. S. Hwang, D. J. Wilson, D. Li, and Y. Jia, “Automated volumetric segmentation of retinal fluid on optical coherence tomography,” Biomedical Optics Express, vol. 7, no. 4, pp. 1577–1589, 2016. [27] H. Bogunović, M. D. Abràmoff, and M. Sonka, “Geodesic graph cut based retinal fluid segmentation in optical coherence tomography,” in Proceedings of ophthalmic medical image analysis, 2015. [28] M. Pilch, K. Stieger, Y. Wenner, M. N. Preising, C. Friedburg, E. M. zu Bexten, and B. Lorenz, “Automated segmentation of pathological cavities in optical coherence tomography scanspathological cavities in oct scans,” Investigative ophthalmology & visual science, vol. 54, no. 6, pp. 4385–4393, 2013. [29] S. Roychowdhury, D. D. Koozekanani, S. Radwan, and K. K. Parhi, “Automated localization of cysts in diabetic macular edema using optical coherence tomography images,” in 2013 35th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC). IEEE, 2013, pp. 1426–1429. [30] A. Lang, A. Carass, E. K. Swingle, O. Al-Louzi, P. Bhargava, S. Saidha, H. S. Ying, P. A. Calabresi, and J. L. Prince, “Automatic segmentation of microcystic macular edema in oct,” Biomedical optics express, vol. 6, no. 1, pp. 155–169, 2015. [31] M. Esmaeili, A. M. Dehnavi, H. Rabbani, and F. Hajizadeh, “Threedimensional segmentation of retinal cysts from spectral-domain optical coherence tomography images by the use of three-dimensional curvelet based k-svd,” Journal of Medical Signals and Sensors, vol. 6, no. 3, pp. 166–171, 2016. [32] G. R. Wilkins, O. M. Houghton, and A. L. Oldenburg, “Automated segmentation of intraretinal cystoid fluid in optical coherence tomography,” IEEE Transactions on Biomedical Engineering, vol. 59, no. 4, pp. 1109– 1114, 2012. [33] W. Wieclawek, “Automatic cysts detection in optical coherence tomography images,” in Mixed Design of Integrated Circuits & Systems (MIXDES), 2015 22nd International Conference. IEEE, 2015, pp. 79– 82. [34] E. K. Swingle, A. Lang, A. Carass, H. S. Ying, P. A. Calabresi, and J. L. Prince, “Microcystic macular edema detection in retina oct images,” in SPIE Medical Imaging. International Society for Optics and Photonics, 2014, pp. 90 380G–90 380G. [35] E. K. Swingle, A. Lang, A. Carass, O. Al-Louzi, S. Saidha, J. L. Prince, and P. A. Calabresi, “Segmentation of microcystic macular edema in cirrus oct scans with an exploratory longitudinal study,” in SPIE Medical Imaging. International Society for Optics and Photonics, 2015, pp. 94 170P–94 170P. [36] L. Zhang, W. Zhu, F. Shi, H. Chen, and X. Chen, “Automated segmentation of intraretinal cystoid macular edema for retinal 3d oct images with macular hole,” in 2015 IEEE 12th International Symposium on Biomedical Imaging (ISBI). IEEE, 2015, pp. 1494–1497. [37] T. Wang, Z. Ji, Q. Sun, Q. Chen, S. Yu, W. Fan, S. Yuan, and Q. Liu, “Label propagation and higher-order constraint-based segmentation of fluid-associated regions in retinal sd-oct images,” Information Sciences, vol. 358, pp. 92–111, 2016. [38] A. González, B. Remeseiro, M. Ortega, M. G. Penedo, and P. Charlón, “Automatic cyst detection in oct retinal images combining region flooding and texture analysis,” in Proceedings of the 26th IEEE International Symposium on Computer-Based Medical Systems. IEEE, 2013, pp. 397–400. [39] A. Rashno, K. K. Parhi, B. Nazari, S. Sadri, H. Rabbani, P. Drayna, and D. D. Koozekanani, “Automated intra-retinal, sub-retinal and subrpe cyst regions segmentation in age-related macular degeneration (amd) subjects,” Investigative Ophthalmology & Visual Science, vol. 58, no. 8, pp. 397–397, 2017.

[40] J. Kohler, A. Rashno, K. K. Parhi, P. Drayna, S. Radwan, and D. D. Koozekanani, “Correlation between initial vision and vision improvement with automatically calculated retinal cyst volume in treated dme after resolution,” Investigative Ophthalmology & Visual Science, vol. 58, no. 8, pp. 953–953, 2017. [41] M. Yang, “Convergence properties of the generalized fuzzy c-means clustering algorithms,” Computers & Mathematics with Applications, vol. 25, no. 12, pp. 3–11, 1993. [42] M. P. Fay and M. A. Proschan, “Wilcoxon-mann-whitney or t-test? on assumptions for hypothesis tests and multiple interpretations of decision rules,” Statistics surveys, vol. 4, p. 1, 2010.

Abdolreza Rashno received the B.S. degree in computer engineering from Shahid Chamran University of Ahvaz, Iran, in 2009 and the M.S. degree (with honors) in artificial intelligence engineering from Kharazmi University of Tehran, Iran, in 2011. Currently, he is a Ph.D. candidate of computer engineering at Isfahan University of Technology, Iran and visiting researcher student at University of Minnesota, USA. His research interests are image segmentation, medical image analysis, image restoration, content-based image retrieval, evolutionary computing and machine learning. Dara D. Koozekanani received the Ph.D. degree in biomedical engineering in 2001, and the M.D. degree in 2003 from Ohio State University, Columbus, OH, USA. His research dissertation involved the application of computer vision techniques to the analysis of optical coherence tomography images. He completed the ophthalmology residency at the University of Wisconsin, Madison, WI, USA, in 2007, and completed a surgical retinal fellowship at the Medical College of Wisconsin, Madison, in 2009. He is currently an Associate Professor of ophthalmology on the clinical faculty at the University of Minnesota, Minneapolis,MN, USA. He sees patients with a variety of surgical and medical retinal diseases. His research interests include the application of ophthalmic imaging technologies and automated analysis of those images. Paul Drayna graduated magna cum laude from Hamline University in St. Paul, MN with a B.A. in biology and history in 2004. He graduated with alpha omega alpha honors from the Uniformed Services University of the Health Sciences F. Edward Hebert School of Medicine in Bethesda, MD in 2008. Dr. Drayna completed his residency training in Ophthalmology at the San Antonio Uniformed Services Health Education Consortium in 2012 and spent the following 3 years as a Comprehensive Ophthalmologist training residents in San Antonio; this included a deployment to Afghanistan to serve as Theater Ophthalmology Consultant. He graduated in 2017 from his fellowship in Vitreoretinal Surgery at the University of Minnesota in Minneapolis, MN. He works as a full-time Vitreoretinal surgeon in the Air Force teaching residents in San Antonio, TX. His research interests include clinical applications of ophthalmic imaging technologies. Behzad Nazari received the B.Sc., M.Sc., and Ph.D. degrees in electrical engineering from Sharif University of Technology, Tehran, Iran, in 1993, 1995, and 2004, respectively. He is currently an Assistant Professor with the Department of Electrical and Computer Engineering, Isfahan University of Technology, Isfahan, Iran. His current research interests include image and signal processing..

Hossein Rabbani received the B.Sc. degree in Electrical Engineering (Communications) from Isfahan University of Technology, Isfahan, Iran, in 2000 with the highest honors, and the M.Sc. and Ph.D. degrees in Bioelectrical Engineering in 2002 and 2008, respectively, from Amirkabir University of Technology (Tehran Polytechnic), Tehran, Iran. In 2007 he was with the Department of Electrical and Computer Engineering, Queen’s University, Kingston, ON, Canada, as a Visiting Researcher, in 2011 with the University of Iowa, IA, United States, as a Postdoctoral Research Scholar, and in 2013-2014 with Duke University Eye Center as a Postdoctoral Fellow. He is now an Associate Professor in Biomedical Engineering Department and Medical Image and Signal Processing Research Center, Isfahan University of Medical Sciences, and Senior Member of IEEE. His main research interests are medical image analysis and modeling, statistical (multidimensional) signal processing, sparse transforms, and image restoration. He has published more than 110 papers and book chapters in these areas.

Keshab K. Parhi (S’85-M’88-SM’91-F’96) received the B.Tech. degree from Indian Institute of Technolgy, Kharagpur, India, in 1982, the M.S.E.E. degree from the University of Pennsylvania, Philadelphia, PA, USA, in 1984, and the Ph.D. degree from the University of California at Berkeley, Berkeley, CA, USA, in 1988. He has been with the University of Minnesota, Minneapolis, MN, USA, since 1988, where he is currently a Distinguished McKnight University Professor and Edgar F. Johnson Professor in the Department of Electrical and Computer Engineering. He has published over 600 papers, is the inventor or co-inventor of 29 patents, has authored the textbook VLSI Digital Signal Processing Systems (New York, NY, USA: Wiley, 1999), and coedited the reference book Digital Signal Processing for Multimedia Systems (Boca Raton, FL, USA: CRC Press, 1999). His research interests include the VLSI architecture design and implementation of signal processing, communications and biomedical systems, error control coders and cryptography architectures, high-speed transceivers, stochastic computing, secure computing, and molecular computing. He is also currently working on intelligent classification of biomedical signals and images, for applications such as seizure prediction and detection, schizophrenia classification, biomarkers for mental disorders, brain connectivity, and diabetic retinopathy screening. Dr. Parhi has served on the Editorial Boards of the IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS PART I AND PART II, the IEEE TRANSACTIONS ON VLSI SYSTEMS, IEEE TRANSACTIONS ON SIGNAL PROCESSING, the IEEE SIGNAL PROCESSING LETTERS, and the IEEE Signal Processing Magazine, and served as the Editor-in-Chief of the IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS PART I from 2004 to 2005. He currently serves on the Editorial Board of the Journal of Signal Processing Systems (Springer). He has served as the Technical Program CoChair of the 1995 IEEE VLSI Signal Processing Workshop and the 1996 Application Specific Systems, Architectures, and Processors conference, and as the General Chair of the 2002 IEEE Workshop on Signal Processing Systems. He was the Distinguished Lecturer of the IEEE Circuits and Systems Society from 1996 to 1998. He served as a Board of Governors Elected Member of the IEEE Circuits and Systems Society from 2005 to 2007. He is the recipient of numerous awards including the 2017 Mac Van Valkenburg award, the 2012 Charles A. Desoer Technical Achievement award and the 1999 Golden Jubilee medal, from the IEEE Circuits and Systems society, the 2013 Distinguished Alumnus Award from IIT Kharagpur, the 2013 Graduate/Professional Teaching Award from the University of Minnesota, the 2004 F. E. Terman award from the American Society of Engineering Education, the 2003 IEEE Kiyo Tomiyasu Technical Field Award, and the 2001 IEEE W. R. G. Baker Prize Paper Award.

Saeed sadri received the B.Sc. and M.Sc. degrees in electrical engineering from Tehran University, Tehran, Iran, in 1976 and 1978, respectively and the Ph.D. degree in electrical engineering from Isfahan University of Technology, Isfahan, Iran, in 1998. He is currently a Professor with the Department of Electrical and Computer Engineering, Isfahan University of Technology, Isfahan, Iran. His current research interests include image and signal processing. View publication stats