Diabetic macular edema grading in retinal images using vector quantization and semi-supervised learning

2.1Overview

The proposed screening system for grading of DME presented can be schematically described by the block diagram in Fig. 1. It consists of three main processing stages: 1) macula localization, containing several main retina tissue detection, i.e. vessel, optic disc and macula; 2) exudate detection, containing three steps of exudate candidates’ identification, feature extraction and suspicious candidates classification; 3) grading of diabetic macular edema, according to macular coordinates and location of exudates to classify the input image into one of the status: normal, Non-CSME and CSME. Once the exudates and macular coordinates are obtained, the grading of DME can be done according to the grading criteria shown in Table 1. The following sections mainly describe the detection methods of macula and exudate in detail.

2.2Macula localization

In the retinal image, macula is the darkest circular region and locates at a specific distance from optic disc (OD). According to the criterion, many methods are proposed [10]. However they fail to detect the macula when the retina contains large dark lesions, i.e. hemorrhages. In this study, a method based on considering multiple features is proposed and achieves a better performance. It is critical to detect the retinal blood vessels and the optic disc prior to localizing the macula.

The method [11] was adopted to assess the size of the OD, and a morphological top-hat operation combined with thresholding technique is employed to segment the main vessels of the retina in the green channel. Then, the OD is detected by using a line operator filter. Following this, a candidate region of interest is defined as the ring area with a radius of 2 to 3 times as the optic disc diameter (DD) from the center of the optic disc [10]. Within the ring area, a circle template with a diameter of about 2/3 of DD [12] is designed and applied to select macula candidates by template matching. Furthermore, two statistic values in the ROI, the mean intensity (M𝑅𝑂𝐼) and the standard deviation of the intensities (σ𝑅𝑂𝐼), are introduced for macula candidate identification. The products of M𝑅𝑂𝐼 and σ𝑅𝑂𝐼 are sorted in ascending order and the ROIs with the top three values are identified as macula candidates. The true macula is selected according to the fact that no blood vessels locate inside the macula. The results of the aforementioned procedure are illustrated in Fig. 2a–e.

Figure 3.

The procedure of the exudate candidates’ identification.

2.3Vector quantization for exudate candidates’ identification

In retinal images, which are usually two-dimensional 24-bit color images obtained through a digital fundus camera, exudates generally appear as bright regions with variable brightness, sizes, locations and shapes. Accurate exudate segmentation is a challenging task due to the large variety in size, intensity, shape and contrast. In this study, the vector quantization (VQ) technique is introduced to segment the exudates. VQ [13] is an efficient technique for data compression, and it can be defined as a mapping function that maps k-dimensional vector space into a finite set CB={C1,C2,C3,…,CN}. The set CB is called a codebook consisting of N code vectors and each code vector Ci={ci⁢1,ci⁢2,ci⁢3,…,ci⁢k} is k-dimensional. It works by dividing a large set of points into groups having approximately the same number of points closest to them if these points are large enough to ensure the probability distribution of points. The groups are normally generated by clustering algorithms, and each group is represented by its centroid as in K-means. In our method, the Linde-Buzo-Gray (LBG) algorithm is used to generate the codebook. More specifically, an initial vector C0 is set as the average of the entire training vectors. By adding a small constant error ϵ>0, the vector is then split into two vectors C10=(1+ϵ)⁢C0 and C20=(1-ϵ)⁢C0. The iterative algorithm is run with the two vectors as the initial codebook, and the process is repeated until the desired size of the codebook is obtained.

Intensity inhomogeneity and uneven illumination hinder the segmentation of the exudates. To solve it, a local region-based segmentation is adopted to segment the image into a number of homogeneous sub-images before the clustering. A region-based segmentation scheme is presented in Fig. 3. Firstly, the OD is removed from the green channel image. Then, the image is partitioned into patches with same size by the grid of M*M. Subsequently, the VQ clustering is carried out on each patch to obtain the codebook and the image is reconstructed using the codebook. The grid of M*M for dividing the image into non overlapping patches is empirically set to 7 * 7, and the size of codebook is also empirically set to 8. Merging the results of all the patches is done followed by selecting brightest region and binarization on each patch. The segmentation process by VQ is defined as Eq. (1).

where I𝑠𝑒𝑔 denotes the exudate result of the whole image by VQ, V⁢(⋅) represents the VQ segmenting to the exudate for the input patch, Pi is the i patch of the all non-overlapping patches and n is the number of all patch. Finally, the bright ROI close to the main vessels are removed by the method in [1], since some bright structures such as the reflections and nerve fibers within the vessels have similar appearance to exudates and can lead to high false positive rate. The final result of exudate candidates is illustrated in Fig. 2f.

2.4Semi-supervised learning for exudate candidates classification

Features extraction is critical for classification. Since exudates appear as yellow or white areas with different size and shape in the retina, color features are extracted from gray level image, the G as well as R channels of RGB color space and the L channel of CIElab color space. Exudates have strong and sharp edges, thus the mean, standard deviation, minimum, maximum of the gradient are calculated on the candidate regions to capture the local gradient information. Entropy is also calculated by including all pixels in bounding square around candidate region to measure the context information of candidate regions. In addition, another two contextual features, the distance to the closest candidate and the number of the nearby candidates, based on the spatial relation with surrounding similar lesions are also taken into account since exudates often occur in clusters and a particular slightly suspect candidate has a higher chance of being a true exudate when another more obvious exudate is present nearby. In total, a 23 dimensional feature is extracted for each candidate exudate.

After extracting features for all candidates, the following phase is to classify them into exudates or non-exudates by machine learning. The traditional classification methods require a sufficiently large number of training datasets with label information. Nevertheless, manually annotating labels for these exudate candidates is difficult. There is no such a publicly available database containing both ground truth of DME diagnosis at an image label and exudate segmentation at a lesion label. Although training a supervised classifier can be done on another datasets with ground truth of exudate, the difference between datasets may lead to a decline in the performance of DME grading. Fortunately, semi-supervised learning (SSL) can provide a suitable solution for this problem. SSL aims to utilize unlabeled data and the limited available labeled data to learn better supervised models with good generalization ability. This study adopt a graph-based semi-supervised learning approach in [9] which defines a similarity graph where vertices are the labeled as well as unlabeled examples, and edges indicate the similarities of pairs of vertices. Let V be a set of vertices with n=|V|, and the edge weights matrix is defined by W∈Rn×Rn. Given two vertices xi and xj, the edge weight Wi,j indicates the similarity of these instances and is calculated using Gaussian RBF kernel: Wi,j=exp⁡(γ⁢||xi-xj||2), where γ is a kernel parameter. With the help of the graph structure, the unlabeled vertices can be accurately classified by using a small amount of labelled examples. Given a labelled instances (xi,yi) and a unlabeled instances xj, where i∈(1,…,L) and j∈(L+1,…,L+U), and a graph 𝑮 is denoted as a (L+U)×(L+U) correlation matrix. The graph-based semi-supervised learning method assumes that a function f can be found so that it is close to the given labels on the labeled examples and it is also smooth on the whole graph. This framework is based on semi-supervised learning with graph embeddings, and a feed-forward neural networks is used to represent embeddings of instances. Given the input feature vector x, the output function of the k-th hidden layer of the network is defined as hk⁢(x)=𝑅𝑒𝐿𝑈⁢(ak⁢hk-1⁢(x)+bk), where ak and bk are parameters of the k-th layer, and h0⁢(x)=x, 𝑅𝑒𝐿𝑈⁢(x)=max⁡(0,x). Given a graph 𝑮, the objective function f of this framework can be written as:

where f𝑆𝐿 is a supervised loss of predicting the labels, f𝑈𝐿 is an unsupervised loss of predicting the graph context, and λ is a regularization parameter.

The first term is the log conditional probability for labeled data, which is defined as:

where w represents the weighting parameter for the class label y. The representation hk⁢(x) from the k layers on the input feature vector x and the representation hl⁢(x) from the l layers on the embedding are concatenated, and fed into a softmax layer to predict the class labels of the instances.

The unsupervised loss is usually defined with the log loss -logp(c|i). Given an instance i and its context c, the objective of the unsupervised loss is formulated as minimizing the log loss of predicting the context c using the embedding of an instance as input features. It uses the embedding of a vertex to predict the context in the graph, where the context is generated by random walk. The distribution p⁢(i⁢c⁢δ) is defined for a sampling precess and is conditioned on labels y1:L and the graph 𝑮. A sampling (i⁢c⁢δ) is done from the distribution. Given (i⁢c⁢δ), the cross entropy loss of classifying the pair (i⁢c) to a binary label δ can be minimized as Π(i,c,δ)⁢log⁡μ⁢(δ⁢aCT⁢ei), where Π⁢(i⁢c⁢δ) is an indicator function that outputs 1 when the argument is true, otherwise 0, μ⁢(⋅) is the sigmoid function, a′⁢s represent the weighting parameters, and ei is the embedding of instance i. The embedding e is used with l1 layers of the network, and the second term in Eq. (2) can be written as

Therefore, according to Eqs (3) and (4), the function f in inductive learning can be written as follow.

In the optimization phase, stochastic gradient descent (SGD) is applied to train the above model and the result images of suspicious exudates classification are shown as Fig. 2g and h.