Automated segmentation of intraretinal cystoid macular edema based on Gaussian mixture model

1. Introduction

Diabetic retinopathy (DR) is one of the main microvascular complications of diabetes, and is the leading cause of vision loss in many developed countries.¹ About half of all patients with DR will develop diabetic macular edema (DME). DME is the most common cause of vision loss among people with DR and it can occur at any stage of the disease course¹ and is characterized by the growth of blood vessels from choroid into the macular region. Hence, early discovery and treatment of DME is imperative. The standard treatment for these diseases is the use of anti-angiogenic intravitreal injections which suppress neovascularization and lead to fluid resorption.² The volume of intraretinal fluid is used as a primary factor for guiding these injections, especially by subjective estimation based on limited spectral domain optical coherence tomography (SD-OCT) slices, which potentially leads to substantial inconsistency in treatment. Furthermore, manual segmentation is time-consuming. Hence, accurate segmentation of the fluid in the retina has an important role in the diagnosis of these diseases.^3,4

OCT, especially SD-OCT, is one of the rapid developing imaging techniques.^3,5,6,7,8 OCT has the advantages of being noninvasive and providing a cross-sectional, micro-scale description of the optical reflectance properties of tissues, which allows the differentiation of retinal structures in the axial direction, enabling clinicians to diagnose retinal diseases better. The development of this technique has revolutionized the noninvasive management of eye diseases.

Given its wide application in ophthalmology, many algorithms based on OCT images have been proposed to segment intraretinal fluid regions including semi-automatic and fully automatic.

Fernandez et al. developed a semi-automatic method to delineate the boundaries of retinal fluid areas using a deformable model, but this method requires the initialization of snake model by human beings.⁹ Zheng et al. developed a semi-automatic algorithm to segment intraretinal and subretinal fluid in OCT which achieved quantitative analysis of the fluid regions. In this approach, an expert clicks on each of the desired candidates after the coarse segmentation and fine segmentation.¹⁰ Wang et al. developed an interactive method to segment the fluid-associated region in 3D SD-OCT slices,¹¹ but this method relies on the initialization by experts.

Since then, many algorithms have been proposed to achieve automatic segmentation. The methods reported in Refs. 12 and 13 are both based on the hard constrain of intensities of whole images. The study reported in Ref. 14 proposed a 3D analysis method to identify the fluid-filled regions, which relies on the segmentation of 10 layers of the retina. The method described in Ref. 15 achieved automated segmentation of fluid regions using neutrosophic sets and graph algorithms, but its key step is the middle-layer segmentation which is easily affected by pathological changes. Generally, it is difficult to segment all retinal layers in OCT images for various diseases. The study presented in Ref. 16 proposed an iterative high-pass filtering approach to localize the cysts, but this method is an extension of layer segmentation. Rashno et al.¹⁷ proposed a method in which a Kernel graph cut is used in neutrosophic domain to segment fluid region. Chiu et al.¹⁸ proposed a method in which segmentation of closed-contour features in ophthalmic images is based on graph theory and dynamic programming. Another study used an automatic method combining region flooding and texture analysis to detect cysts.¹⁹

Due to the recent developments in machine learning, more supervised and unsupervised methods have been proposed to segment retinal fluid regions.

Many automatic segmentation methods based on supervised learning have been developed to segment the fluid regions, such as k-nearest neighbor classifier,^20,21,22 kernel regression,²³ random forest classifier,^24,25 and Adaboost classifier.²⁶ A new three-dimensional curvelet transform-based dictionary learning for automatic segmentation of intraretinal cysts was also proposed.²⁷ But the performance of the supervised methods^{20,21,22,23,24,25,26,27} relies on the selection of labels. Several methods based on deep learning have also been proposed. A study²⁸ proposed a convolutional neural network to segment macular edema. Venhuizen et al.²⁹ used two FCNNs to segment the inner retinal cysts (IRC). A deep learning encode–decode model was also built to segment the intraretinal cystoid fluid.³⁰ Roy et al.³¹ built a new framework named ReLayNet to segment the retinal layer and fluid. Tennakoon et al.³² proposed an encoding-decoding framework inspired by U-net to segment retinal fluid. Girish et al.³³ proposed an FCNN method for vendor-independent IRC segmentation. Gopinath et al.³⁴ proposed a method based on FCNN by selective enhancement to segment retinal cysts. Bogunović et al.³⁵ used a machine learning approach combined with an effective segmentation framework based on geodesic graph cut to segment retinal fluid regions. Although the methods described in Refs. 28–35 achieved good performance, the algorithms used require devices with high hardware performance.

Unsupervised methods have also been proposed for segmentation of retinal fluid. Wang et al. proposed an automated volumetric segmentation method to detect retinal fluid on OCT image using fuzzy level-set,³⁶ but its performance on images with unclear boundaries and irregular shape of edema was not sufficient.

For OCT images with irregular shape edema and unclear boundaries (shown in Fig. 2(b)), coarse segmentation and refinement of boundaries are useful. In this paper, we propose a novel automatic algorithm based on Gaussian mixture model (GMM) and level-set. First, the proposed method utilizes the GMM clustering to obtain the coarse mask and then uses the adjustable thickness-distance information to get an accurate mask for initialization of level-set. To the best of our knowledge, this is the first study to use GMM to segment fluid regions on OCT images. Finally, the region-based active model is employed to refine the segmentation. Spatial information is also introduced to search for the potential fluid. Figure 1 shows the flowchart of our method. The proposed method is validated on the challenging dataset using DME from patients with DR, which demonstrates that the proposed method has a good performance on OCT images with irregular shape edema and unclear boundaries (shown in Fig. 2(b)).

2. Methods

2.1. Data acquisition

In this paper, we use data comprising 10 SD-OCT cubes from 10 eyes of 10 patients diagnosed with DR, resulting in a total of 1280 B-scans. All cubes were obtained using a Cirrus SD-OCT cube from 128 contiguous $512\times 1024$ pixel B-Scan images (each B-Scan comprising 512 A-Scans containing 1024 pixels).

2.2. Pre-processing

SD-OCT imaging is based on the interferometric detection of coherence optical beams which have low temporal coherence and high spatial coherence. Thus, SD-OCT images are obtained with speckle noise. Speckle in SD-OCT tomograms is influenced by several factors, such as the wavelength of the imaging beam and the structural details of the imaged object.³⁷ The speckle sizes differ in the axial and lateral dimensions, which are mainly determined by the source bandwidth and numerical aperture, respectively. In addition to speckle noise, shot noise also exists in SD-OCT images, which is additive in nature and can be adequately described by the additive white Gaussian noise process.³⁸ To solve this problem, we apply a modified bilateral filtering algorithm³⁹ to reduce noise. Compared with the conventional denoising methods, such as the patch-based methods^40,41 or the Gaussian filtering, this method is better for edge preservation and has a relatively low time complexity. Moreover, layer segmentation is an important auxiliary step for retinal cystoid segmentation, because the cysts are located on the retinal region. Therefore, a restricted region is defined within the inner limiting membrane (ILM) and Bruch’s membrane for segmentation of cysts. A publicly available three-dimensional segmentation software, Iowa reference algorithm^42,43 is used to obtain two layers: ILM layer and inner boundary of retinal pigment epithelium (IB-RPE) layer. Results of the layer segmentation on normal retina and DR cases are shown in Fig. 2, indicating that ILM layer and IB-RPE layer are segmented efficiently.

2.3. Cysts segmentation

First, the intensity values of each pixel in each B-scan are clustered by GMM clustering to obtain coarse segmentation.⁴⁴ Second, the thickness-distance threshold is utilized to remove false positive regions. Finally, the optimized level-set method is applied to refine the boundaries and spatial information is utilized to search for potential fluid regions.

It should be noticed that the assumption of Gaussianity and stationarity for the signal are important in many aspects of statistical modeling in the field of signal processing. The standard methods in signal processing are justified based on the theory of Gaussian Stationary Processes (GSP).⁴⁵ As reported in Ref. 46, the distribution of OCT intensities in the retina should be a mixture of exponentials. In Ref. 47, Gaussianization was used to convert the probability distribution function (pdf) of each OCT intra-retinal layer to a Gaussian distribution. Hence, it is reasonable to model the pdf of an image such as OCT image with Gaussian function. The GMM clustering is a soft clustering method widely used in many fields with good results including image segmentation and moving object detection in video sequences.⁴⁸ In this paper, we postulate that the distribution of OCT intensities is the mixture of $k$ Gaussian functions ($k=2$). This hypothesis means that the whole image can be divided into two clusters, one is region of interest and the other is the remaining regions. As shown in Fig. 3(b), the two green lines show the intensities distribution of regions of interest and whole B-scan intensities. Hence, the intensity distribution of each B-scan can be seen as the sum of $k$ components including target region and other regions. The function can be defined as follows :

where ${\alpha }_i$ is the probability of $i$th Gaussian function, $x$ represents the intensity of each pixel, ${\mu }_i$ and ${\sigma }_i$ are the parameters of $i$th Gaussian function. In OCT images, cysts are similar to vitreous fluid and choroid regions due to their dark or low intensity. The intensity probability distribution of OCT image is shown in Fig. 3, and the distribution of the whole image is regarded as a mixture of two Gaussian distributions (shown as Fig. 3(b)). GMM clustering algorithm groups all pixels into $k$ different clusters based on the posterior probability of pixel belonging to each cluster. In this work, we divided the pixels in each OCT image into two groups on the basis of the intensity of each pixel. The initialization of these parameters (${\mu }_i$, and, ${\sigma }_i$) depends on the maximal and minimum value of the intensities.

After GMM clustering, the cluster with minimum intensity $C_{min}$ is chosen as the original segmentation for subsequent processing. Given the intensity similarity, retinal fluid vitreous and other low intensity regions in or below retina are included in the cluster $C_{min}$. When the macular edema is seen in the retina, the areas over ILM layer and below the IB-RPE layer are removed. The cluster classification map which is processed by region restriction is shown in Fig. 3(c). The retinal fluid vitreous and other low intensity regions below retina are removed based on the layer segmentation of ILM and IB-RPE. The region of interest (ROI) (Fig. 3(c)) is located between ILM and IB-RPE layers.

2.4. Restriction of the target region using adjustable threshold and thickness information

The outer nuclear layer (ONL) contains the nucleus of cone cell and rod cell. Because of the high liquid content in the nucleus, the signal is low in OCT image, which limits effective segmentation (shown in Fig. 3(c)). Figure 3(c) shows that the ONL presented in the map has a banding distribution along the IB-RPE layer. The distribution of distance between each pixel and IB-RPE layer is shown in Fig. 4, and there is a peak distribution between 0 and 100, which is identical to the distance distribution of ONL band. This distribution is used to restrict the target region. Given that the distribution is similar to Gaussian function, this is leveraged to threshold the distance: (1) Calculate the mean value of the first peak. (2) Flip-duplicate the left part of the curve using the mean value. (3) Find the zero point of the new curve. The distance value at the zero point is chosen as the threshold. Then, the false lesion region is removed from the coarse segmentation, as shown in Fig. 4.

There are segmentation artifacts (shown as the red batch in Fig. 5) similar to the real retinal fluid in size, intensity and shape. To improve the accuracy of the algorithm, we use thickness information to remove the false segmentation. As the fluid builds up, the thickness of retina gradually enlarges. The thickness information described in the function is defined as follows :

where $d$ is the distance of each pixel to the IB-RPE layer, thick is the thickness between ILM layer and IB-RPE layer of each batch in the map, Mean($d)$ is the average distance. The diagrammatic drawing is shown in Fig. 5. For this step, we set threshold $r<4$ and Mean($d)<M$ for each pixel. If the pixel satisfies these two conditions simultaneously, the pixel will not be regarded as fluid. The image without thickness thresholding and the resulting image are shown in Fig. 6.

2.5. Refinement of the boundaries

Level-set methods are widely used in the image segmentation process, including nature image segmentation and medical image segmentation. In recent years, this method has been applied to detect and segment abnormality in medical images. Given that the quality of our data is not good, the segmentation did not fit the boundaries well. In this work, a region-based active contour model⁴⁹ is applied to refine the boundary of retinal fluid. The segmentation results after removing artifacts are used as the initialization of this level-set method.

2.6. Search for potential fluid

To overcome the break between frames, the spatial information is used to search for potential fluid. According to OCT imaging theory, adjacent B-scans show adjacent parts of retina. Furthermore, the edemas have similar locations in adjacent B-scans. Thus, the connection within consecutive five slices is used to recognize the missing fluid areas, defined as

where mask ($i$) is the $i$th B-scans from 128 B-scans. According to the different weights of adjacent images, the score of each pixel belonging to the fluid can be calculated. The nearer the adjacent slice is from target image, the higher its weight is. If the score is no less than 7, the pixel is regarded as fluid. The adjacent images do not need to be fully aligned because Eq. (3) is used to search the intersection of consecutive 5 slices. The outputs of this step are used as the initialization seeds for the level-set method. Furthermore, the level-set method will evolve the accurate boundaries and will delete the false regions according to the seeds. An example is shown in Fig. 7.

3. Results

In this paper, 10 cubes are used to assess the accuracy and robustness of the proposed method. We compared automatic segmentation results with ground truth contoured by experienced graders to show the correlation between the results obtained by the proposed method and those of experts. A comparison was also made between our method and the fuzzy level-set algorithm³⁶ to evaluate the performance of the proposed method in images with edema with unclear boundaries. As the proposed method is based on unsupervised clustering, we also added the comparison with classical unsupervised clustering — $K$-means clustering combined with level-set and the FCM–level-set method.³⁶ To evaluate the contribution of the adjustable distance-thickness threshold combined with spatial information, we compared the proposed method with using GMM clustering and region restriction (GMMR) approaches. The performance of the proposed method is measured from qualitative and quantitative experiments.

3.1. Qualitative results

Wang’s method³⁶ is used to test our datasets and to compare the performance of this method and the proposed method. In Wang’s method,³⁶ they used images acquired by Optovue SD-OCT system, of which the image is a real three-dimensional with high quality. The sweep spacing of Cirrus device is longer than that of Optovue due to differences in the SD-OCT system. Hence, the voting of three directions is not suitable to our data, and the FCM-level-set segmentation in one direction is therefore used for comparison. At the same time, the GMMR method is used to test our datasets to compare the benefits of our method with Wang’s method. Figure 8 shows the comparison among the three methods and manual segmentation. Figure 9 shows a three-dimensional visualization of fluid segmentation using the proposed method.

3.2. Quantitative results

Three criteria are utilized to evaluate the segmentations quantitatively, including dice similarity coefficient, true positive error and false positive error. The Dice similarity coefficient is as follows :

where $V_a$ stands for our automatic segmentation results, $V_m$ is the ground truth that is contoured by experienced grader. The Dice coefficient ranges from 0 to 1, where 0 suggests the two results are completely different and 1 denotes that they are similar. True positive error defined as Eq. (5) denotes the rate of correctly detected area. False positive error defined as function Eq. (6) indicates the fraction of incorrectly detected area.

Quantitative experiments are conducted to determine how the proposed method compares with the manual grading. In Table 1, the mean and standard deviation of TPVF, FPVF and DSC for segmentation by GMM, GMMR, GMMR–level-set and our proposed method are determined to evaluate the effect of each step of the proposed method. GMM compared with the proposed method means that we only use the GMM method to achieve fluid segmentation. GMMR indicates that firstly the GMM is used to achieve the coarse segmentation, and then the areas above ILM layer and below IB-RPE layer are removed. GMMR–level-set stands for the combination of GMMR and level-set. Because the gray scale distribution of the background is similar to that of the edema, many background pixels are mistaken as edema, resulting in very low accuracy, so we introduce the area limitation and thickness constraint to improve the accuracy. By comparing GMM and GMMR, we can summarize that the use of region constraint helps to improve the DSC and depress the FPVF. The use of level-set method greatly improved the TPVF, which could be known by comparing GMMR and GMMR–level-set. Although the mean TPVF of the proposed method is slightly lower than that of GMM, the propose method is generally superior to the other three algorithms, so this conclusion could be summarized by its lowest FPVF, highest DSC. Therefore, we introduced the region limitation and thickness constraint to improve the accuracy. The value for DSC in GMM (Table 1) is small because the results by GMM include the areas above the ILM layer and below the IB-RPE. The FCM–level-set method and $K$-means–level-set method are also compared with our proposed method shown as Fig. 10. The $K$-means–level-set method is the combination of $K$-means and level-set, and contains the step of removing areas which are above the ILM layer and below IB-RPE layer. From left to right, the three pictures represent the box plots of DSC, TPVF and FPVF, respectively. According to Fig. 10, the proposed method shows better stability than the other three methods. Due to the distance and thickness restrictions, real fluid regions may be omitted by mistake in the proposed method, hence the mean TPVF of the proposed method is lower compared with that of the FCM–level -set, but our method has the smallest standard deviation. The mean DSC of our method is higher than the other three algorithms and it has the smallest standard deviation. Although the TPVF in $K$-means–level-set is much higher, the DSC in $K$-means–level-set is the lowest. As a result, the proposed method shows good performance and stability in intraretinal cystoid macular edema segmentation. Due to the low quality of our data, the DSC of the three methods is low. Figure 11 shows the consistency between manual segmentation and our method. Analysis of the Bland–Altman plots reveals that the 100% limits of agreement were [$-0.22$, 0.07] for the proposed algorithm and that of the expert.

4. Discussion

This paper proposes a novel method to automatically segment DME in low quality images, which includes four steps: (1) Denoising the image and segmenting retinal layers; (2) Recognizing retinal fluid using GMM clustering; (3) Restricting the region using thickness-distance threshold; (4) Refining the segmentation using level-set and spatial information. The GMM clustering step is used to classify all pixels in each B-scan so that the candidate pixels are clustered. The novelties can be summarized as follows: (1) It is the first attempt to use GMM clustering to segment retinal edema. The results show that it is effective in segmenting OCT images. (2) We introduce the adjustable threshold to constrain the ROI independent of the hard constrain. (3) Spatial information between consecutive five frames is used to search for the potential fluid. These steps are combined to realize the whole segmentation. In this method, multiple parameters are applied including intensity, thickness-distance threshold and spatial information which overcomes the limitations of using image intensities only, and the use of morphological operation. The distance information is used by fitting the distance distribution.

The experimental results demonstrate that the proposed method can detect the intraretinal cystoid macular in low quality images and is comparable with manual grading. Figures 8, 2(c), 2(d), 4(c) and 4(d) show that the FCM–level-set method and GMMR mistook the artifacts as the fluid. Although all the three methods achieved the segmentation of retinal fluid, the boundaries of FCM method shown in Fig. 8(2b) and GMMR shown in Fig. 8(4b) are not smooth. The proposed method matches smoothly with the manual segmentation. The quantitative results in Table 1 show the performance and the effect of each step of the proposed method. Figure 10 shows the performance and stability of the proposed method compared with FCM–level-set, $K$-means–level-set and GMMR. Figure 11 shows the consistency between the proposed method and manual segmentation.

One limitation of the proposed method is that it relies on the layer segmentation. Although we used ILM and IB-RPE layers, the accuracy of these two layers may influence the results. Segmentation of IB-RPE may affect the threshold of distance. As the disease progresses, some colobomas may occur in ILM which will obscure the boundary of ILM, thus the fluid below the ILM will be regarded as vitreous fluid. This will increase the false negative rates and decrease the DSC.

Our group plans to improve the GMM method to achieve more accurate segmentation results. We shall do this by using more features like texture and location information to improve the specificity of fluid areas. In addition, we would leverage the advantages of features at different resolution. Finally, we plan to refine the iteration function to recuse the cost of time, at the same time enhance the accuracy of initialization.

Conflict of Interest

The authors declare that there are no conflicts of interest related to this article.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (NSFC) (Grant Nos. 61701192, 61471226 and 61671242), the Natural Science Foundation of Shandong Province, China (Youth Fund Project) under Grant No. ZR2017QF004, Natural Science Foundation of Shandong Province (Grant Nos. JQ201516 and 2018GGX101018), the Taishan scholar project of Shandong Province (No. tsqn2016023), Fundamental Research Funds for the Central University (30920140111004), and the China Postdoctoral Science Foundation (No. 2017M612178).