Fig. 1
An example of an original RGB fundus image of the healthy eye and individual colour channels of the image. In standard fundus image, the red (R) channel appears oversaturated, while the green (G), and the blue (B) channels show the blood vessels and retinal nerve fiber layer highly contrasted
Fig. 2
An example of OCT volume and circular scans. a SLO image (left) with the volume scan pattern allocated by the green lines and one B-scan (right) measured at the position depicted by the blue line in SLO image; b SLO image (left) with the circular scan pattern defined by the blue circle and the B-scan (right) measured along this circle in direction given by the arrow. The curves in individual B-scans define segmentation of the RNFL
3 Methods
An illustrative schematic diagram of the proposed RNFL assessment methodology is depicted in Fig. 3. The texture analysis is carried out within the peripapillary area at the locations without the blood vessels only. Our previously published matched filtering approach [19] is used for the blood vessel segmentation. Various regression models are tested towards prediction of the RNFL thickness using the proposed texture features. The regression models are trained on small square image regions (ROIs) selected from fundus images in the database and known measurement of the RNFL thickness. Circular profiles are extracted from the predicted images provided by the regression models. The resulted profiles are further validated with respect to the real RNFL thickness measured via OCT. The following subsections deal with the description of particular processing steps as well as evaluation of the approach.
Fig. 3
Schematic diagram of the proposed methodology
3.1 Data Preprocessing
3.1.1 Preprocessing of Fundus Images
The fundus images are preprocessed in several steps. First, standard uncompressed TIFF format is reconstructed from the raw data, whereas a linear gamma transfer function is applied in the reconstruction process. Secondly, non-uniform illumination of fundus images is corrected together with the increase of image contrast using CLAHE (Contrast Limited Adaptive Histogram Equalization) technique [20]. The RNFL texture is the most contrasted in the green (G) and the blue (B) channels of the input RGB image (Fig. 1). Therefore, an average of G and B channel (called GB image) is computed for each fundus image after CLAHE. Further, only the GB images are used for processing.
In the first step, we manually selected small square-shaped image regions of interest (ROIs) with size of 61 × 61 pixels from all fundus images included in the group of normal subjects. Extraction of ROIs was performed uniformly in the peripapillary area to the maximum distance not exceeding 1.5 × diameter of the OD; whereas only locations without the blood vessels were considered (Fig. 4). In this way, a total number of 354 ROIs was collected. Particular ROIs then represent the typical RNFL pattern depending on the position in the peripapillary area for normal subjects without any signs of glaucoma disease.
Fig. 4
At the top: section of the GB image after CLAHE processing with depiction of ROI positions; at the bottom: few examples of magnified ROIs with the RNFL texture taken in different positions in the peripapillary area (around the OD)
3.1.2 Preprocessing of OCT Data
The OCT volume data were preprocessed in order to obtain the RNFL thickness in the peripapillary area of each subject in the database. Hence, the RNFL was segmented and the corresponding RNFL thickness map was created using freely available research software [21]. Segmentation of the RNFL layer was done automatically with very good precision so that only subtle manual corrections had to be performed in some B-scans using this software package (see segmentation of the RNFL in Fig. 2), especially in the area of large blood vessels (shadow artifacts in the B-scans). The final RNFL thickness map can be seen in Fig. 5.
Fig. 5
The RNFL thickness map mapped on the SLO image of a normal subject. The colour spectral scale represents the changes of RNFL thickness approx. from 20 µm (red) to 180 µm (green)
3.1.3 Fundus-OCT Image Registration
Our previously published [22] landmark-based retinal image registration approach with manually selected landmarks and second-order polynomial transformation model was applied for registration of fundus to OCT–SLO image data. This registration step was necessary to be able to compare the proposed texture features with the RNFL thickness at various positions on the retina. However, different registration approaches could be used for this purpose as well, e.g. as in [23].
3.2 Feature Extraction
The advance texture analysis methods, namely Gaussian Markov random field (GMRF) [24] and local binary patterns (LBP) [25] were used for the description of RNFL texture. These approaches were selected due to their rotation- and illumination-invariant properties as well as noise robustness.
3.2.1 Gaussian Markov Random Fields
First set of features is given by GMRF non/causal two/dimensional autoregressive model [24]. The model assumes the image texture is represented by a set of zero mean observations y(s) [24]:
(1)
for a rectangular M × M image lattice Ω. An individual observation is then represented by the following difference Eq. [24]:
(2)
where N s is a neighborhood set centered at pixel s, ɸr is the model parameter of a particular neighbor r, and e(s) is a stationary Gaussian noise process with zero mean and unknown variance σ. A neighborhood structure depends directly on the order and type of the model. We assume a fifth-order symmetric rotation-invariant neighborhood structure as shown in Fig. 6. The structure considers five parameters expressed by particular numbers. These five parameters describe the relationship between central pixel and its neighbors. Gaussian variance σ is the sixth parameter of the model. Then, these 6 parameters represent features, which are used for the RNFL texture description.
Fig. 6
A fifth-order symmetric rotation-invariant neighborhood structure
The least square error (LSE) estimation method is used for estimation of the GMRF model’s parameters according to the following equations [24]:
![$$j ={{\left[ \sum\limits_{\Omega }{q(s){{q}^{T}}(s)} \right]}^{-1}}\left( \sum\limits_{\Omega }{q(s)y(s)} \right),$$](https://i0.wp.com/basicmedicalkey.com/wp-content/uploads/2017/06/A329170_1_En_2_Chapter_Equ3.gif?w=960)
(3)
(4)
where
![$$q(s)=col\left[ \sum\limits_{r\in {{N}_{i}}}{y(s+r)};\ i=1,...,I \right],$$](https://i0.wp.com/basicmedicalkey.com/wp-content/uploads/2017/06/A329170_1_En_2_Chapter_Equ5.gif?w=960)
(5)
for an i-th–order neighborhood structure.
3.2.2 Local Binary Patterns
The second applied method—LBP is based on conversion of the local image texture into the binary code using rotation-invariant and uniform LBP operator [25]. The local image texture around the central pixel (x c , y c ) can be characterized by the LBP code, which is derived via the Eq. [25]:
(6)
where U(G p ) means:
(7)
In Eqs. 6 and 7, g c corresponds to the grey value of the central pixel (x c, y c) of a local neighborhood and g p (p = 0,…, P–1) corresponds to the grey values of P equally spaced pixels on a circle of radius R (R > 0) that form a circularly symmetric neighborhood structure. The LPB operator expressed by Eg. 6 assumes uniform patterns. The “uniformity” of a pattern is ensured by the term U(G P ). Patterns with U value of less than or equal to two are considered as “uniform” [25]. It means these patterns have at most two 0–1 or 1–0 transitions in the circular binary code.
Two variants of LBP were utilized in the proposed approach. Both variants are based on the rotation-invariant and uniform LBP16,2 operator (i.e. P = 16, R = 2). One variant uses only LBP distribution computed from an input GB image. Then, the grey-level histogram of such parametric image is computed and extraction of 6 statistical features follows [25]: mean value, standard deviation, skewness, kurtosis, total energy and entropy. In the second variant, standard LBP distribution is supplemented with computation of local contrast C P,R :
(8)
Then, in turn, a joint histogram of 
and 
(LBP/C) is computed. A feature vector is then obtained from LBP/C joint histogram by extraction of 14 texture features proposed by Haralick et al. [26] and Othmen et al. [27] (energy, contrast, homogeneity, entropy, correlation, sum average, sum variance, sum entropy, difference variance, difference entropy, two information measures of correlation, cluster shade, and cluster prominence).
and (LBP/C) is computed. A feature vector is then obtained from LBP/C joint histogram by extraction of 14 texture features proposed by Haralick et al. [26] and Othmen et al. [27] (energy, contrast, homogeneity, entropy, correlation, sum average, sum variance, sum entropy, difference variance, difference entropy, two information measures of correlation, cluster shade, and cluster prominence).3.2.3 Pyramidal Decomposition
Finally, a 26-dimensional feature vector assembled via connection of particular texture analysis approaches (GMRF, LPB, and LBP/C) is obtained. In addition, the features are computed for an original image resolution and even for each of the two levels of Gaussian pyramid decomposed images. Let the original image be denoted as G 0(i,j), which is zero level of the Gaussian pyramid. Then, the l-th level of the pyramid is defined as follows:
(9)
where w(m,n) is a two-dimensional weighting function, usually called as “generating kernel”. According to, [28] recommended symmetric 5 × 5 kernel, written in separated form as ![$$w=\left[ \frac{1}{4}-\frac{a}{2},\frac{1}{4},a,\frac{1}{4},\frac{1}{4}-\frac{a}{2} \right]$$](https://i0.wp.com/basicmedicalkey.com/wp-content/uploads/2017/06/A329170_1_En_2_Chapter_IEq3.gif?w=960)
, where a = 0.4, is utilized. Finally, a 78-dimensional feature vector is obtained via extraction of the features from G 0(i,j), G 1(i,j), and G 2(i,j). Composition of the final feature vector is depicted schematically in Fig. 7.
, where a = 0.4, is utilized. Finally, a 78-dimensional feature vector is obtained via extraction of the features from G 0(i,j), G 1(i,j), and G 2(i,j). Composition of the final feature vector is depicted schematically in Fig. 7.
Fig. 7
Schematic diagram of the final feature vector
3.2.4 Feature Selection and Regression
The aim of this work is to propose the utilization of texture analysis in fundus images for description of changes in the RNFL pattern related to variations in the RNFL thickness. The ability of the proposed texture analysis methods, in connection with several regression models, to predict the RNFL thickness has been investigated. Different regression models—linear regression (LinReg) [29], two types of support vector regression (v-SVR, ε-SVR) [30], and multilayer neural network (NN) [31] have been tested to predict values of the RNFL thickness using the proposed texture features. In addition, different feature selection approaches [32] have been tested towards identification of the most relevant feature subset of the original feature set. Finally, we have chosen a popular wrapper-based feature selection strategy with a sequential forward search method (SFS) that provided the most accurate prediction of the RNFL thickness using various regression models. In SFS strategy, standard forward hill-climbing procedure is utilized. The procedure starts with an empty feature set and sequentially adds a feature that yields in the best improvement of a subset. This proceeds until there is no improvement in performance of a particular feature subset. In each iteration of the wrapper approach (for each feature subset), a cross-validation procedure is used to evaluate model output via a chosen evaluation criterion (e.g. mean squared error).
4 Results and Discussion
4.1 Evaluation Methodology
Spearman’s rank correlation coefficient (ρ) and root mean squared error of prediction (RMSEP) were used as evaluation criteria of the models output. ρ is computed between the model predicted output y and the target variable c as follows [33]:
