2.2.1 Preprocessing.

The purpose of preprocessing and contrast enhancement is to remove any artifacts that occur during retinal image acquisition process. We applied the following procedures to preprocess the raw images.

1. Convert the colored (RGB) image to green channel image and process the contrast and enhancement. Perhaps one of the most basic enhancement techniques is the contrast stretch. We applied the decorrelation stretching method for the enhancement works by increasing the differences in hue [14], this technique gives better results for microvasculature in retinal image (provide more vessels' information).
2. We then applied the hybrid median filtering method to separate vessels from others. This approach is the windowed filter of nonlinear class that can easily remove impulse noise while preserving edges (during the vessel segmentation). We first generates ranking in multiple sub-neighborhoods for each pixel and then takes the median values from those rankings and performs a final ranking to select the final result. This procedure better preserves fine details while removing the noise. More details can be seen in [15].

In image analysis, Non Uniform Illumination in an image often leads to diminished structures and inhomogeneous intensities due to different texture of the object surface and shadows cast from different light source directions. It is one of the most important factors affecting the appearance of an image, more adverse in case of biological images. For the correction of non uniform illumination, we apply the morphological operations. The basic idea is as follow: display the background approximation as a surface first, and then under the assumption that the vessels maintain predominantly in local linear orientation we select a line structure to generate a morphological filter for vessel segmentation [16]. This approach was used to solve the problem by estimating background approximation as a surface in order to extract the non-uniform background from the image, and then construct the new image by subtracting this estimated background from the original image. Results reveal that various particles present in the image with exact boundaries along with the removal of non-uniform illumination at the background. Basic morphological operators are erosion, dilation, opening, closing as defined in [16]. The top-hat, bottom-hat, and morphological enhancement operators are written as(1)Because blood vessels are darker than the background in retinal image, a bottom-hat operator and enhancement operator are used to generate the feature response images. After this, we applied the dual-tree complex wavelet (CWT) on the above enhanced image since it is a valuable enhancement of the traditional real wavelet transform that is nearly shift invariant and, in higher dimensions (is oriented in 2D here), directionally selective [17]. At this point, the image denoising task was completed and we have eliminated most of the small or rough artifacts that may not be detected and removed by the above two approaches. This is especially beneficial for the vessels extraction. Other alternative way is to apply the matched filters approach on this enhanced image. An example of the above preprocessed result is presented in Fig. 4 (sample images from a public database of MESSIDOR).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 4. Retinal image after preprocessing; (a) original green channel image; (b) enhanced/preprocessed image; (c) vessels extracted image.

https://doi.org/10.1371/journal.pone.0075699.g004

In the next section we will discuss features extraction and formulation. There are three approaches for the texture analysis: 1) Statistical texture analysis involving gray level co-occurrence matrix (GLCM) and run-length matrix, 2) HOS analysis and 3) fractal analysis to generate useful features from the enhanced image. There are a total of 75 features generated initially. It includes 50 features coming from HOS analysis and 17 features generated from statistical texture analysis, and 8 features (three from mono-spectrum and five from multi-spectrum) from fractal analysis.

2.2.2 Features extraction and formulation - Statistical based texture analysis.

All related texture features from the statistical-based texture analysis were generated from each retinal image that derived from gray level co-occurrence matrix (GLCM) and run-length matrix [18]–[26]. For an image of size M×N, the gray level co-occurrence matrix (GLCM) is defined as in [20](2)Where (p, q), (p+Δx, q+Δy)∈M×N, d = (Δx, Δy) and |⋅| denotes the cardinality of a set I(p,q) indicates the gray level of a pixel at (p,q) in an image. Thus given a gray level i in an image, the probability that the gray level of a pixel at a (Δx, Δy) distance away is j is(3)From each co-occurrence matrix we computed the following features:Energy measures the textural uniformity of the image, i.e., the repetition of pixel pairs. It is the measurement of the denseness or order in the image. Entropy measures disorder or randomness of the image and it is an indication of the complexity within an image, thus, more complex images have higher entropy values. Contrast is a measure of the presence of local variations (or differences in the GLCM) in the image, and higher contrast values indicate large local variations. Homogeneity (also called an inverse difference moment) is inversely proportional to the contrast at constant energy. Similarly at constant contrast, homogeneity is inversely proportional to energy [27]. It measures how close the distribution of elements in the GLCM is to the diagonal of GLCM.

Other measurements such as Moments 1–4 are defined as:(4)where g is the integer power exponent that defines the moment order. Moments are the statistical expectation of certain power functions of a random variable and are characterized as follows [20]: moment 1 is the mean which is the average of pixel values in an image [21]; moment 2 is the standard deviation; moment 3 measures the degree of asymmetry in the distribution; and moment 4 measures the relative peakedness or flatness of a distribution and is also known as kurtosis [22].

Other similar approach can also be used: having the GLCM normalized, we can then derive eight second order statistic features which are also known as haralick features [23] for each image, which are: contrast, correlation, energy, entropy, homogeneity, dissimilarity, inverse difference momentum, maximum probability. In addition to these features, we also applied correlation, dissimilarity, inverse difference momentum and maximum probability, which is different from above mentioned features.

The gray level run-length matrix (RLM) is defined as the numbers of runs with pixels of gray level i and run length j for a given direction [24]. RLMs was generated for each sample image segment having directions (0°,45°,90° & 135°), then the following five statistical features were derived: short run emphasis, long run emphasis, gray level non-uniformity, run length non-uniformity and run percentage [26]. Basically it allows extraction of higher order statistical texture features.

2.2.3 Features extraction and formulation - High order Spectra analysis (HOS).

Higher order spectra (HOS) are known to have the ability to detect non-linearity and deviations from Gaussian assumption. Motivated by these, a set of HOS based parameters were proposed as features to differentiate the new vessels from normal vessels. HOS are spectral representations of higher moments and they are derived from the averaged Fourier spectrum signal. The bispectrum B(f₁, f₂), of a signal is the Fourier transform (FT) of the third order correlation of the signal. It is given by(5)where X( f) is the FT of the signal x(nT), * represents complex conjugation and E[.] stands for the expectation operation. It retains Fourier phase information. The frequency f may be normalized by the Nyquist frequency to be between 0 and 1. The bispectrum, given by equation (5), is a complex-valued function of two frequencies. The bispectrum which is the product of three Fourier coefficients exhibits symmetry and need only be computed in a non-redundant region. Assuming that there is no bispectral aliasing, the bispectrum of a real-valued signal is uniquely defined with the triangle 0≤f₂≤f₁≤f₁+f₂≤1.This is termed as Ω, the principal domain or the non-redundant region (See the triangle region in Fig. 5).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 5. Non-redundant region (Ω) of computation of the bispectrum for real signals. Parameters are calculated from this region.

https://doi.org/10.1371/journal.pone.0075699.g005

Briefly, HOS based and spectral based features are: Mean of spectral magnitude for HOS [28]:(6)where B( f₁, f₂) is the bispectrum of the signal.

Other features such as Entropy 1, 2 and 3 are defined as:And the feature of bispectrum phase entropy (EntPh):whereφ: Bispectrum phase angle

L: Number of points within the samples

We used each of these bispectral invariant features for every 18° from 0° to 180°. Therefore we obtained total 50 features of HOS. Next we will discuss the model-based methods such as fractal model approach. The fractal model is useful for modeling certain natural textures that have a statistical quality of roughness at different scales, and also for texture analysis and discrimination.

2.2.4 Features extraction and formulation - Fractal analysis (FA).

Fractal model (analysis) may be considered as the model-based texture analysis. We can classify it into mono-fractal and multi-fractal. Fractals are geometric objects whose increasing details under magnification resemble exactly or statistically the whole object (self-similarity). Such fractal objects are not easily “measurable” in classic geometric terms because some of their physical characteristics (length, mass, area, volume, and so on) are largely dependent on the magnification used when they are measured. This means that the surface of the fractal object is often complex and any patterns presented on the surface may be combined with different physical characteristics. The diversified complexity of fractal may be described with the concept of fractal dimension [29], which may easily describe the incompleteness or fragmentation of an entirety. Moreover, recent studies show that such surface complexity of image may be described not only by its fractal dimension but also its multifractal spectra, which is mathematical description of a surface that can accurately reflect its features, and is compatible with the various theoretical models that related to surface structures. Thus we can apply fractal analysis to determine the surface complexity which is measured on the gray scale, and by using multifractal spectra one can obtain more detailed information than is possible with the fractal dimension alone. Recently this kind of technique is widely used in retinal vessels analysis.

Recalling fractals are characterized by scale invariance and it always shows a similar degree of covering space (shape). The way in which scale invariance appears is characteristic for the considered structure and can be expressed by a single value, the fractal dimension (FD). FD is a non-integer value which reflects the structure's “convoluteness”. An excellent introduction of fractal analysis can be found elsewhere [30]–[31]. It has also been speculated that fractal analysis could be helpful in diagnosing diabetic retinopathy [32].

Since the retinal vasculature is a fractal that follows the theory of fractal geometry [30], [32]–[33], the fractal dimension is expected to be a natural measure of new vessel formation. It has also been shown that the new vessel formation changes the fractal dimension of the vessel pattern. The fractal dimension exhibits a high degree of sensitivity with respect to new vessel formation, while the image preprocessing with respect to the representation of the individual vessel thickness did not affect the results. The fractal dimension appears to be the “natural” measure for proliferative changes and could be used for automated detection of proliferative diabetic retinopathy in the future [34].

Measuring fractal dimension has previously been attempted to quantify small changes to the human retinal vasculature, not immediately apparent by human observation, and act as an early marker of disease [35]–[37]. Mono-fractal analysis is an indicator of vascular change, it has achieved limited success as retinal vessels may have different characteristics found on different location or with different scale of measurement. Greater success has been reported by considering the retinal vascular pattern to be multi-fractal, characterized by a hierarchy of exponents rather than a single fractal dimension [38].

We first applied box-counting algorithm (with Hausdorff dimension) approach to calculate fractal dimension in binary type of image (segmentation of vessels) and then the Fourier fractal dimension (FFD) approach, which have been proposed by [35]–[37]. FFD has been used to quantify the grayscale images projected on to 3-D fractal surface [38]. The advantage of FFD is that it computes the fractal dimension of gray scale images, and eliminates the need for image segmentation [29]. It has also been found to be relatively insensitive to noise and it is believed to work effectively with data having low signal-to-noise ratio [38]–[39]. We adapted the similar FFD approach proposed by Azemin M. and et al. [37]. The parameters such as slope and intercept were generated based on this approach. Finally we applied the multifractal spectra since it can describe the evolution of the probability distribution of fractal structures [35], [38]–[40]. Instead of using simple fractals (or monofractals), multifractals are characterized by a hierarchy of exponents, rather than a single fractal dimension. In this paper the box-counting method was used to characterize multiple spectra. [38].

Multifractals main parameter is Hólder's exponent [41]–[42].Where μ(box) represents dimension of box and ε dimension of longitude of box.

Notes: Local and global information from the spectrum are used for segmentation, noise deletion and edge detection at picture points. Segmentation is an important step for description of the basic individual process. One of the commonly used approaches is filtering, but this approach has a main disadvantage on the lost of precision due to preliminary filtering.

Alternative approach is observing of the image as measure of the fix resolution. Irregularities of this measure can be examined with the help of multifractal analysis. The general principle is as follow:

1. - First, different measures and capacities are defined from image with a gray level.
2. - After that corresponding multifractal spectrum has to be calculated, enabling local (using α) and global (using f(α)) information. No hypothesis about the signal regularity was used.

Importance and advantage of fractal and multifractal analysis (MFA), compared to “classical” signal analysis is one way of handling of irregularities. MFA tries to extract the information directly from the singularity, whereby “classical” approach often observes LF (low-pass filter) filtered versions, with different filtering depths for irregularities observing and noise repressing. Based on specific values of α and f(α), a homogeneity point can be isolated in original signal. By image pixel extraction, which satisfy chosen value of the parameter α or spectrum f(α), it is possible to extract by any of the known methods. Additional advantage is that such segmentation causes no degradation of the starting image: all pixels interrelations stay unchanged and therefore image details are being completely kept. In order to describe fast changing signals (small area) of retinal image and to express variability, we need to examine fractal characteristic. Usage of classic statistical methods in such case (mean values) can cause valuation errors. Significant singularities are indicated by multifractality of the process.

It is also concluded that the fractal dimension is a descriptor of early changes of pathological vessels changes, which introduce further understanding of the mechanisms of complex changes. Fractals are very important in both medicine and pathology. In our study, five parameters were generated from multifractal analysis: three from local (α) and two from global (f(α)).

Based on the above three approaches, the features set formulated from the texture analysis is given as I_t = {τ₁, τ₂, …, τ_n}; where n = 75 is the number of features set for each image.

2.3.1 Dimension reduction.

We first reduce the dimension since some of the features may be correlated or redundant in the model. A penalized logistic regression method (Pelora) [43] was used to cluster the datasets from the above procedures combining all 75 parameters (factors). This is a supervised clustering algorithm that has been used with external information about response variables for clustering genes in genetic study. This algorithm is mainly based on penalize logistic regression analysis and it combines feature selection, supervision, feature clustering and sample classification in a single step. This approach has an advantage of dealing with unbalanced problem (outcomes) since penalized likelihood is a general approach to reducing small-sample bias in maximum likelihood estimation [44]. On the other hands, this approach is different from the classical logistic regression in that it groups features into clusters instead of extracting individual features. This method was used for features selection with the purpose of dimension reduction generated from retinal images and then classified the variables (potential factors) obtained. This approach will preserve the property of the classified (identified) groups with complex interactions.

2.3.2 Model selection with multimodel inference.

After the dimensions have been reduced, we applied the method of automated model selection and model-averaging that provides a wrapper for GLM and similar functions, automatically generated all possible models with the specified response of exudates(s) and explanatory variables, and determine the best models with a defined criterion (e.g. AIC or AICc). The best model is mainly based on the bias-variance trade-off.

Statistical models with probability density functions were used, where we can maximize the likelihoodby minimize the Kullback-Leibler distance from the true model to the approximating model, i.e.,Here KL distance is one way to measure the distance between two densities:(7)g( y) for the discrepancy is weighted by the probability observing the data y, and g( y)/f( y,θ) is the ratio that measures the discrepancy between the two models for data y. AIC selects the model closest to the true model on average. Since the true distribution of g is unknown, so we replace it with the empirical (observed) distribution. Thus the estimation of KL distance is introduced and a bias will also be given. To determine the expected value of this bias, the Takeuchi Information Criterion (TIC) was used. Finally under the assumption of correct approximate model, we derived AIC and related criteria.

In order to deal with the model selection uncertainty and usual sampling uncertainty, we applied AIC to select the model closest to the true model on average (Atypical datasets have lower influence than typical ones). To do this, we first apply bootstrap technique to resample the data based on our raw data and then produce the distribution of best models. For instance, after boostrap data we apply the smoothed model weights and finally determine: 1) the value of the parameter on average (also its variance) for solving the model selection uncertainty problem and, 2) the standard error for a given model that follows a distribution (conditional variance of the estimator), this is to deal with usual sampling uncertainty. Thus we can combine the above two variances to form the overall variance for the estimator . Such approach reduces the spurious estimates (false positive) problem, which standard model selection techniques often encountered (see Freedman's paradox [45]).

We then applied GLMULTI [46], an R package for automated model selection and multi-model inference with GLM and related functions. The basic idea of this approach is that from a list of explanatory variables, GLMULTI builds all possible models involving these variables and, optionally, included their pair-wise interactions (for the computing of simplicity and easy interpretation reasons we will not consider higher order interaction in this paper). Restrictions can be specified for candidate models, by excluding specific terms, enforcing marginality, or controlling model complexity. Detailed information is referred to [46],[47]–[49].