An ensemble of classifiers based on different texture descriptors for texture classification

2.1 Local quinary patterns (LQP)

The LQP operator represents an improvement of the classical LBP operator (Ojala et al., 2002). The main weakness of the LBP operator is that the binary $$s(x)=\left\{\begin{array}{c}1\text{,}x⩾0\hfill \\ 0\text{,}x<0\hfill \end{array}\right.$$

function thresholds exactly at the intensity value of the central pixel q_c, which makes this operator sensitive to noise especially in near-uniform regions. To overcome this weakness, the quinary coding proposed in Paci et al. (2012) can be used. The quinary coding is achieved introducing two thresholds τ₁ and τ₂ in the canonical LBP s(x) function which becomes: $$S(x\text{,}{\tau }_1\text{,}{\tau }_2)=\left\{\begin{array}{c}\hfill 2\text{,}x⩾{\tau }_2\hfill \\ \hfill 1\text{,}{\tau }_1⩽x<{\tau }_2\hfill \\ \hfill 0\text{,}-{\tau }_1⩽x<{\tau }_1\hfill \\ \hfill -1\text{,}-{\tau }_2⩽x<-{\tau }_1\hfill \\ \hfill -2\text{,}x<-{\tau }_2\hfill \end{array}\right.$$

To reduce the size of the histogram summarizing the extracted patterns, the original LQP code is thus split into 4 different sets of binary patterns, according to the following binary function b_c(x), c ∈ {−2,−1, 1,2}: $$b_c(x)=\left\{\begin{array}{c}\hfill 1\text{,}x=c\hfill \\ \hfill 0\text{,}\mathit{otherwise}\hfill \end{array}\right.$$

Each set of binary patterns is mapped according to the “uniform rotation invariant” (riu2) mapping Ojala et al., 2002 using two neighborhoods: (i) the 8 pixel neighborhood resulting in 10 features and (ii) the 16 pixel neighborhood resulting in 18 features. So the histogram obtained by combining the two neighborhoods is composed by 28 elements.

Finally, these 4 histograms are concatenated to form the feature vector containing 112 bins: 28 (features of each histogram) × 4 (number of histograms).

2.2 Local phase quantization (LPQ) and LPQ with ternary coding (LPQ3)

LPQ is based on the blur invariance of the Fourier Transform Phase (Rahtu et al., 2012), locally computed from the 2D Short Term Fourier Transform (STFT) for each pixel position of the image over a square neighborhood.

Considering the 2D Fourier transforms F of the original picture f, H of the point spread function (PSF) h and G of the blurred image g, they are bounded by the following relation $$G=F\cdots H$$

Thus magnitude and phase aspects can be separated thus resulting in $$|G|=|F|·|H|$$

and $$\angle G=\angle F+\angle H$$

If the blur PSF h is centrally symmetric, its Fourier transform H is always real-valued, and its phase is a two-valued function given by $$\angle H=\left\{\begin{array}{cc}0\text{,}\hfill & H⩾0\pi \text{,}H<0\hfill \end{array}\right.$$

meaning that F and G have the same phase at those frequencies which make H positive.The LPQ method is based on the above properties of blur invariance. It uses the local phase information extracted using STFT computed over a square neighborhood at each pixel position x of the image f(x). The binary code assignment at each x depends on the phase information only and the phase is computed by observing the component of the F_x vector $${\mathbf{F}}_x=[\text{Re}\{{\mathbf{F}}_x^c\}\text{,Im}\{{\mathbf{F}}_x^c\}]\text{,}$$ where ${\mathbf{F}}_x^c=[F({\mathbf{u}}_1\text{,}\mathbf{x})\text{,}F({\mathbf{u}}_2\text{,}\mathbf{x})\text{,}F({\mathbf{u}}_3\text{,}\mathbf{x})\text{,}F({\mathbf{u}}_4\text{,}\mathbf{x})]$ and the u vectors frequency vectors are u₁ = [a,0]^T, u₂ = [0,a]^T, u₃ = [a, a]^T, and u₄ = [a,−a]^T where a is small enough to satisfy H(u) ⩾ 0.We defined G_x = V^TF_xwhere V is an orthonormal matrix derived from the singular value decomposition of the covariance matrix of the transform coefficient vector F_x. The resulting vectors are quantized: $$q_j=\left\{\begin{array}{c}\hfill 1\text{,}{g}_j⩾0\hfill \\ \hfill 0\text{,}{g}_j<0\hfill \end{array}\right.$$ where g_j represents the j-th component of G_x. These quantized coefficients are represented as integers between 0 and 255 using the binary coding $$b={\displaystyle \sum _{j=1}^8}{q}_j{2}^{j-1}$$

These integer values are then organized in the feature vector useful for classification tasks. In this paper we have used a rectangular neighborhood of size 3 and 5 concatenating the two histograms before training a given classifier.

As for LBP a non-binary coding is proposed for the LPQ operator (Paci et al., 2012): in this case the LPQ with ternary coding (LPQ3) uses the following quantizer: $$q_j=\left\{\begin{array}{cc}\hfill 1\text{,}\hfill & \hfill g_j⩾\rho ·{\tau }_j\hfill \\ \hfill 0\text{,}\hfill & \hfill -\rho ·{\tau }_j⩽g_j⩽\rho ·{\tau }_j\hfill \\ \hfill -1\text{,}\hfill & \hfill g_j⩽-\rho ·{\tau }_j\hfill \end{array}\right.$$ where τ_j is set to half of the standard deviation of the j-th component of G_x, and ρ is the given weight. The quantized coefficients are then represented as integers in the interval 0–255 using the following binary codings $$b_{+}={\displaystyle \sum _{j=1}^8}(q_j==1)·2^{j-1}$$

and $$b_{-}={\displaystyle \sum _{j=1}^8}(q_j==-1)·2^{j-1}$$

b₊ and b₋ values are then summarized in two distinct 256 bins histograms and the two histograms are then concatenated thus providing a 512 valued feature vector (for both the neighborhood of size 3 and 5, i.e. the final feature vector is 1024 bins), useful for classification tasks.

2.3 The multi-threshold approach

The use of ternary and quinary codings requires respectively one (τ) or two (τ₁ and τ₂) thresholds to be set. The threshold selection is a critical task in order to reduce the sensitivity to noise of these new operators: thresholds were usually set manually in order to get the best performance in specific problems, but some automatic adaptive procedures were also proposed in Vécsei et al. (2011) exploiting local statistics as mean value and standard deviation inside each neighborhood. Another approach lies in choosing a set of different thresholds for the ternary coding (or a set of different couples of thresholds for the quinary coding), in order to (i) extract a different feature set according to each threshold (or couple of thresholds), (ii) use each feature set to train a different classifier and (iii) fuse all these results according to a fusion rule (sum, mean, vote, etc.).

In details, in this work we considered:

• 25 threshold couples (τ₁ = {1,3,5,7,9} and τ₂ = {τ₁ + 2, …, 11}) for LQP i.e. 25 feature sets;

• 5 thresholds (τ ∈ {0.2, 0.4, 0.6, 0.8, 1}) for LPQ, i.e. 5 feature sets.

This led to new texture descriptors such as the Multi-threshold Local Quinary Patterns (MLQP) and the Multi-threshold Local Phase Quantization with ternary coding (MLPQ3) Paci et al., 2012.

Since we gathered 30 feature sets we trained 30 different SVMs using one feature set each. For the testing part, after computing the same 30 feature sets from the testing images, each SVM classifies the testing feature set corresponding to the same feature set it was trained with. Then the 30 partial classification results are fused together according to the sum rule (Bianconi et al., 2007). Let us define d_t,j(x) the similarity of the pattern x to the class j obtained using the classifier t. The values of d_t,j(x) are: real values that represent the distances from the margin when SVM is the classifier; real values varying in the range 0;1 when GPC is the classifier (see Section 3). Before the fusion the set of similarities of each classifier is normalized to mean 0 and standard deviation 1 (Kuncheva, 2004). The sum rule among a pool of T classifiers is the sum among the T set of scores i.e. Sum rule: ${\mu }_{\text{j}}(\mathbf{x})={\sum }_{t=1}^T{\text{d}}_{\text{t,j}}(\mathbf{x})$ . Each membership to a class is a sum of real values, so it is almost impossible that a given pattern has the same membership to two (or more) classes¹. Each pattern x is assigned to the class with higher membership, i.e. the class of x is argmax_j μ_j(x).

3.1 Support vector machines

Support Vector Machines (SVMs) were first introduced in Vapnik (1995) and are maximum margin classifiers. They are two-class classifiers that find a decision surface by projecting the training samples in the multidimensional feature space and by maximizing the distance of the closest points in the training set to the decision boundary. The goal is to build the hyperplane that maximizes the minimum distance between two classes. Under the hypothesis of linearly separable data $$\mathbf{w}·\mathbf{x}+b⩾1\text{,}{\text{x}}_{\text{i}}\in \text{Class}\_1\text{;}$$ $$\mathbf{w}·\mathbf{x}+b⩽-1\text{,}{\text{x}}_{\text{i}}\in \text{Class}\_2\text{;}$$

the hypothesis space is summarized in H: fw, b = sgn(w⋅x + b).

In order to maximize the distance $$d({\mathbf{x}}_{\mathbf{i}}\text{;}\mathbf{w}\text{,}b)=\frac{|{\mathbf{x}}_{\mathbf{i}}·\mathbf{w}+b|}{\Vert \mathbf{w}\Vert }$$

between the training samples and the hyperplane, the function $$\mathrm{\Phi }(\mathbf{w})=\frac{1}{2}\Vert \mathbf{w}{\Vert }^2$$

should be minimized under the constraints $$y_i(\mathbf{w}·\mathbf{x}+b)⩾1$$

The final form of the decision hyperplane is: $$f(\mathbf{x})={\displaystyle \sum _{i=1}^k}{\alpha }_i{y}_i{\mathbf{x}}_i·\mathbf{x}+b=\mathbf{w}·\mathbf{x}+b$$ where α_i and b are the solutions of a quadratic programming problem.

Unknown test data x_t can be classified by simply computing: $$f(\mathbf{x})=\mathit{sign}({\mathbf{w}}_0·{\mathbf{x}}_{\mathbf{t}}+b_0)$$

It can be seen by examining the above equation that the hyperplane is determined by all the training data, x_i, that have the corresponding attributes of α_i > 0.

In order to get the minimal number of errors during the classification task the constraints can be relaxed, by using the tolerance parameters ξ_i and the penalty parameter C, in the form $$y_i(\mathbf{w}·{\mathbf{x}}_{\mathbf{i}}+b)⩾1-{\xi }_i$$

thus changing the function Φ to be minimized as follows: $$\mathrm{\Phi }(\mathbf{w})=\frac{1}{2}\Vert \mathbf{w}{\Vert }^2+C\left({\displaystyle \sum _{i=1}^l}{\xi }_i\right)$$

The final form of the classifier does not change while now 0 ⩽ α_i ⩽ C. For more information about SVM, the reader is referred to (Cristianini and Shawe-Taylor, 2000).

We have tested radial basis function and linear kernel, the parameters setting is performed in each dataset. To select the parameters an internal 10-fold cross validation is performed using the training data (so the test set is blind) i.e. each training set is randomly partitioned into 10 equal size subsamples, a single set is retained for testing the model, and the remaining 9 sets are used as a new training set. The cross-validation process is then repeated 10 times (the folds), and the results from the 10 folds are averaged. Using this protocol we select the best parameters for SVM (using a grid search approach as suggested in LibSVM (xxxx)), then these parameters (selected without using the test data) are used to classify test patterns.

3.2 Gaussian process classifiers

A Gaussian process is a generalization of the Gaussian probability distribution. A Gaussian process classifier (GPC) is a discriminative approach where the class membership probability is the Bernoulli distribution (Rasmussen and Williams, 2006). GPC has a proven track record classifying many different tasks.

GPC is based on methods of approximate inference for classification, since the exact inference in Gaussian process models with likelihood functions tailored to classification is intractable. In the rest of this section we use the same symbols and notation from Rasmussen and Williams (2006).

In the basic case of a binary GPC classifier, the main idea consists in introducing a Gaussian Process over the latent function f(x), which is used as the argument of an activation function to get a final output varying in the range 0;1. A commonly used activation function is the logistic function $$\sigma (f(\mathbf{x}))=\frac{1}{1+e^{-f(\mathbf{x})}}$$

thus getting a prior on $$\pi (\mathbf{x})=p(y=+1))=\sigma (f(\mathbf{x}))\text{.}$$

The inference problem is made of two steps:

• computing the distribution of f(x) corresponding to a test case $$p(f_{\ast }|x\text{,}\mathbf{y}\text{,}{\mathbf{x}}_{\ast })=\int \mathbf{p}({\mathbf{f}}_{\ast }|\mathbf{X}\text{,}{\mathbf{x}}_{\ast }\text{,}\mathbf{f})\mathbf{p}(\mathbf{f}|\mathbf{X}\text{,}\mathbf{y})\mathbf{df}$$ where $$p(\mathbf{f}|X\text{,}\mathbf{y})=\frac{p(\mathbf{y}|\mathbf{f})p(\mathbf{f}|X)}{p(\mathbf{y}|X)}$$ represents the posterior over the latent variables.

• this distribution over the latent variables is used to produce the prediction for the test case $${\pi }_{\ast }=p(y_{\ast }=+1|X\text{,}\mathbf{y}\text{,}{\mathbf{x}}_{\ast })=\int \sigma ({\mathbf{f}}_{\ast })=\int \sigma ({\mathbf{f}}_{\ast })\mathbf{p}(\mathbf{fast}|\mathbf{X}\text{,}\mathbf{y}\text{,}{\mathbf{x}}_{\ast }){\mathbf{df}}_{\ast }$$

Due to the non gaussian nature of the posterior p(f|X, y), the first integral is not analitically tractable and approximations have to be used.

The Laplace approximation allows the replacement of p(f|X, y) with its Gaussian approximation $$q(\mathbf{f}|X\text{,}\mathbf{y})=N(\mathbf{f}|\hat{\mathbf{f}}\text{,}{A}^{-1})$$ around the maximum of the posterior $$\hat{\mathbf{f}}={\text{arg max}}_{\mathbf{f}}p(\mathbf{f}|X\text{,}\mathbf{y})$$ where $$A=-\nabla \nabla \log p(\mathbf{f}\Vert X\text{,}\mathbf{y}){|}_{\mathbf{f}=\hat{\mathbf{f}}}$$ Following an iterative process to define $\hat{\mathbf{f}}$ and A and computing the mean and variance for f_∗, the prediction can be written as $${\pi }_{\ast }\approx \int \sigma (f_a\mathit{st})q(f_a\mathit{st}|X\text{,}\mathbf{y}\text{,}{\mathbf{x}}_a\mathit{st}){\mathit{df}}_{\ast }$$ where q() is gaussian.

Due to execution time we have used default parameters in all the tested datasets². For more mathematical details of this classification method, the reader would be best served to refer to chapter 3 in Rasmussen and Williams (2006), electronically available at <http://www.gaussianprocess.org/gpml/>. In the experiments we used the code available at <http://www.gaussianprocess.org/gpml/code/matlab/doc/>.