*Article* **Wiener–Granger Causality Theory Supported by a Genetic Algorithm to Characterize Natural Scenery**

**César Benavides-Álvarez 1,\*,†, Juan Villegas-Cortez 2,†, Graciela Román-Alonso 1,† and Carlos Avilés-Cruz 2,†**


Received: 9 May 2019; Accepted: 21 June 2019; Published: 26 June 2019

**Abstract:** Image recognition and classification have been widely used for research in computer vision systems. This paper aims to implement a new strategy called Wiener-Granger Causality theory for classifying natural scenery images. This strategy is based on self-content images extracted using a Content-Based Image Retrieval (CBIR) methodology (to obtain different texture features); later, a Genetic Algorithm (GA) is implemented to select the most relevant natural elements from the images which share similar causality patterns. The proposed method is comprised of a sequential feature extraction stage, a time series conformation task, a causality estimation phase, causality feature selection throughout the GA implementation (using the classification process into the fitness function). A classification stage was implemented and 700 images of natural scenery were used for validating the results. Tested in the distribution system implementation, the technical efficiency of the developed system is 100% and 96% for resubstitution and cross-validation methodologies, respectively. This proposal could help with recognizing natural scenarios in the navigation of an autonomous car or possibly a drone, being an important element in the safety of autonomous vehicles navigation.

**Keywords:** classification; content–based image retrieval; genetic algorithms; image retrieval; image classification; Wiener-Granger causality

### **1. Introduction**

One of the challenges researchers face today is developing an artificial authentication system that has acquisition and processing capabilities similar to those possessed by humans [1]. Artificial vision is defined as the capacity of a machine to see the world that surrounds it in a 3-Dimensional form starting from a group of 2-Dimensional images [2]. Since there is no effective algorithm that can fully recognize any object one can imagine in the entire environment, computer vision is considered an open problem. A computer vision system is composed of different stages that work together for solving a particular problem [3].

Automatic image recognition is among the problems that might be solved using computer vision systems. Researchers are eager to develop these systems and different techniques have been implemented for their improvement, such as machine learning, pattern recognition and evolutionary algorithms.

One of the tasks of an automated image recognition system is to successfully classify and identify natural scenery images (It is said that a scene is natural if the image has no intervention or alteration by

human hands). Currently, thousands of images are generated via different kinds of sources on a daily basis and the constant increase of the Internet has influenced human life.

More than half of the information on the Internet is images, 85% of which were taken with mobile devices with a final estimation of 5 trillion images reported so far [4].

In order to use this information efficiently, an image recovery system based on Content-Based Image Retrieval (CBIR) is necessary. It will help users to find relevant images based on their self-content features or those which are "seen" to e related to them, from our visual perception, even when there is no previous knowledge of the database, such as manual labeling of the images.

Our previous work successfully applied the CBIR technique to the face recognition problem [5,6]. The multiple textures, objects in unknown positions and their different compositions in natural scenery images challenge the proposals that combine different techniques for obtaining a better performance of natural scenery image classification. In this work, we use CBIR feature extraction as an input of a texture causality engine to characterize 5 scenery types, manually defining a base dictionary conformed by 4 textures. In future work, conforming this dictionary is planned to be dynamical, considering more base textures and scenery types to improve classification performance.

In this work, an image retrieval system of natural scenery images is developed by applying the Wiener-Granger Causality (WGC) theory [7] as a tool for analyzing images throughout self-content information. The causal relationships between local textures contained in an image were identified, leading to characterization of a descriptive pattern of a set of scenes inside an image dataset. The selection of causality relationships was carried out using genetic algorithm (GA) implementation as an evolutionary process.

The major stages involved in the developed system are the following (See Figure 1):


**Figure 1.** Proposed general methodology applied for image recognition.

The paper proposes a causality analysis of the natural scenery classes based on a pre-established texture dictionary and the WGC analysis from the CBIR methodology [5,8] in order to provide a whole dataset characterization.

This approach aims to improve the optimization process of evolutionary algorithms. In this case, since the GA [9] shows a simple and fast implementation, it was employed to select the relationships of the local-texture statistical features handled as time series.

Finally, an improvement in the classification accuracy obtained by our proposed strategy is reported, getting 100% on re-substitution and up to 96% for cross-validation methodologies. This approach was implemented using the computer power of a 19-processor cluster and the MPI parallel programming tool.

The current methodology was probed with two databases of natural scenery:


Visualizing the future implementation of an autonomous system of recognition of natural scenes mounted on a car—which will be managed by our proposal as the autonomous system [12–14]—recognizing natural scenarios in the navigation of an autonomous car or possibly a drone, with a 100% certainty, this proposed system will be an important element in the safety of autonomous vehicles.

The rest of the paper is organized as follows: Section 2 presents the state of the art of the problem of image analysis from the CBIR criterion and the WGC theory used in our project, as well as the theoretical support of the WGC model to be applied; in Section 3, the proposed methodology for applying the WGC theory in the natural scenery image characterization is presented; in Section 4 our GA implementation approach to optimize the selection of texture causality relationships is explained; the parallel implementation of our proposal to get good efficiency when processing a large number of images is provided in Section 5; finally, the results and conclusions are presented in Sections 6 and 8, respectively.

### **2. State of The Art**

The problem of image classification and recognition has been studied with different approaches for supporting visual search for different purposes.

Several techniques have been applied successfully to the face recognition problem [6,15–18]. The solutions are favored by controlling the way in which the images are obtained by determining the amount of light, the orientation, the distance, and so forth, in order to obtain ideal face images. In addition, the points to be identified on a face image are well known. The multiple textures, objects in unknown positions and their different compositions make it quite difficult to recognize and identify natural scenery in an image or group of images.

One of the most recent solutions for the classification of natural scenery is the use of the deep learning technique [19], which consists of a set of neural networks connected with each other in successive layers, where each layer network performs a convolution operation on the information of the previous layer, as we can see in Reference [20]. This methodology has the disadvantage of requiring high-end computational resources (memory and CPU) for the training task, unlike the CBIR technique which can be implemented in systems with few resources.

When using CBIR for scenery image classification, significant descriptors are determined considering the image self texture attributes to have an important and effective recovery. In this system, a user presents an image query and the system returns similar images from the database. In Figure 2, the general diagram of a CBIR-based classification system of natural scenery images is shown.

**Figure 2.** Classical methodology of image classification.

One of the first papers that uses the CBIR methodology for natural scenery classification is that by J. Vogel [21]; this work defines a regular 10 × 10 grid on the image; from each grid coordinate an analysis window is opened. Local information is extracted from a window texture and compared with a base texture dictionary; then, the author defines a classification system for natural scenery. A point to be improved is the definition of the base texture dictionary that is set manually, including only typical textures perceived by the researchers. This approach obtains up to 75% average for the cross validation classification test.

Unlike a grid, in Reference [8], random points are thrown on the image and around each point a window is opened; from each window, statistical texture local information is extracted to be grouped, conforming dynamically to a base texture dictionary. The testing is performed considering the generated dictionary obtaining an 85% classification average of natural scenery data bases.

In Reference [22], the CBIR approach is presented to classify natural scenery images through the composition of relevant features in relation to the texture, like in Reference [23], the shape and distribution of the luminosity.

CBIR, being an unsupervised learning technique, still has some disadvantages, since the information extracted is only treated as a histogram that represents the composition of textures in a scenery. This way of characterizing scenery has not been able to obtain more than 85% classification, that is why new proposals that use hybrid methodologies to give CBIR greater robustness arise [24]

In References [24,25], the authors combine the CBIR information with certain semantic content introducing high-level concept objects, trying to link content-based images to objects extracted inside them. This work obtained a percentage of natural scenery classification not greater than that of References [21] and [8].

In this work, a hybrid method of three components is presented. The basic component is CBIR, which generates the information regarding the local texture features of the image. Unlike performing only statistical management of the obtained features (using histograms), in this proposal the second component is responsible for applying a causality technique based on the Wiener Granger causality theory to identify the causal relationships that exist within the basic textures of a type of scenery. Since the causality

component generates different configurations of causal relationships, the third component consists of a GA that allows the selection of the configuration that obtains the best classification percentage for each scenery.

Evolutionary Computational Vision (ECV) as a research area is currently growing in artificial intelligence through two areas of work—computational vision and evolutionary computation. Beginning from a practical point of view, ECV seeks to design the software and hardware solutions necessary to solve hard computer vision problems [1]. Bio-inspired computation within computational vision contains a set of techniques that are frequently applied to hard optimization problems. Its chief objective is to generate solutions formulated in a synthetic way and the artificial evolutionary process based on the evolutionary theory developed by Charles Darwin is the one frequently applied in Reference [9].

#### *2.1. Theoretical Fundamentals Of Wgc*

The causal inference paradigm has been used in different fields of science, for example, in neurology the WGC theory [26] is used to examine areas of the brain and the causal relationships among them. WGC analysis was carried out using sensors [27,28], and, lately in MRI images [29–31], the WGC theory is being used for the study of causal relationships among areas of the brain. Other science fields where WGC theory has been applied is video processing for indexing and retrieval [32]. Video processing for massive people and vehicle identification [33–35] and complex scenery analysis [36]. In this proposal, for the first time, WGC theory is applied to a natural elements and natural scenes retrieval.

In this section, the theoretical framework of the WGC is established. For simplicity and in order to avoid extending mathematically, the theory is presented only for three random processes, being extendable to *n*−processes. In our approach, a random process corresponds to a signal reading associated to one type of texture within a natural scenery; so, for the present analysis, each texture reading corresponds to one stochastic process represented by *Ti*, being *i* the *i*-th texture which has a stochastic behavior disposed into a scenery.

#### *2.2. Stochastic Autoregressive Model*

We assume that each texture can be represented by an autoregressive model into time series. In the current analysis, we will only carry out with three signals, {*T*1, *T*2, and *T*3}, being easily extendable to *n* signals/textures. Let *T*1, *T*2, and *T*3 be three stochastic processes, individually and jointly stationary. Each stationary process can be represented by an autoregressive model in the following way:

$$T1(t) = \sum\_{k=1}^{\infty} \mathbb{C}\_{T1}^1(k) T1(t - k) + \eta\_{T1'}^1 \quad \text{with} \quad \Sigma\_{T1}^1 = \text{var}(\eta\_{T1}^1), \tag{1}$$

$$T\mathbf{2}(t) = \sum\_{k=1}^{\infty} \mathbb{C}\_{T2}^1(k) T\mathbf{2}(t-k) + \eta\_{T2}^1 \quad \text{with} \quad \Sigma\_{T2}^1 = \text{var}(\eta\_{T2}^1), \tag{2}$$

$$T\Im(t) = \sum\_{k=1}^{\infty} \mathcal{C}\_{T3}^1(k) T\Im(t-k) + \eta\_{T3\prime}^1 \quad \text{with} \quad \Sigma\_{T3}^1 = \text{var}(\eta\_{T3}^1), \tag{3}$$

being *η*<sup>1</sup> *<sup>T</sup>*1, *<sup>η</sup>*<sup>1</sup> *<sup>T</sup>*<sup>2</sup> and *<sup>η</sup>*<sup>1</sup> *<sup>T</sup>*<sup>3</sup> random Gaussian noise with zero mean and unit standard deviation; *<sup>C</sup>*<sup>1</sup> *<sup>T</sup>*1(*k*), *<sup>C</sup>*<sup>1</sup> *<sup>T</sup>*2(*k*) and *C*<sup>1</sup> *<sup>T</sup>*3(*k*) are the coefficients of the regression model for textures *<sup>T</sup>*1, *<sup>T</sup>*2 and *<sup>T</sup>*3, respectively.

The joint autoregressive model for the three textures is defined by the equations:

$$T1(t) = \sum\_{k=1}^{\infty} \mathcal{C}\_{T1}^{1,1}(k)T1(t-k) + \sum\_{k=1}^{\infty} \mathcal{C}\_{T2}^{1,2}(k)T2(t-k) + \sum\_{k=1}^{\infty} \mathcal{C}\_{T3}^{1,3}(k)T3(t-k) + \eta\_{T1\prime}^2$$

$$\text{with } \sum\_{l=1}^{2} = \text{var}(\eta\_{T1}^2) \quad \text{(4)}$$

$$T2(t) = \sum\_{k=1}^{\infty} \mathbb{C}\_{T1}^{2,1}(k)T1(t-k) + \sum\_{k=1}^{\infty} \mathbb{C}\_{T2}^{2,2}(k)T2(t-k) + \sum\_{k=1}^{\infty} \mathbb{C}\_{T3}^{2,3}(k)T3(t-k) + \eta\_{T2}^2$$

$$\text{with} \\ \sum\_{l=2}^{2} = \text{var}(\eta\_{T2}^2) \quad (5)$$

$$T\Im(t) = \sum\_{k=1}^{\infty} \mathcal{C}\_{T1}^{3,1}(k)T\mathbf{1}(t-k) + \sum\_{k=1}^{\infty} \mathcal{C}\_{T2}^{3,2}(k)T\mathbf{2}(t-k) + \sum\_{k=1}^{\infty} \mathcal{C}\_{T3}^{3,3}(k)T\Im(t-k) + \eta\_{T3}^2$$

$$\text{with } \Sigma\_{T3}^2 = \text{var}(\eta\_{T3}^2) \tag{6}$$

where ∑<sup>2</sup> *<sup>T</sup>*1, Σ<sup>2</sup> *<sup>T</sup>*<sup>2</sup> and <sup>Σ</sup><sup>2</sup> *<sup>T</sup>*<sup>3</sup> are the variance of the residual terms *<sup>η</sup>*<sup>2</sup> *<sup>T</sup>*1, *<sup>η</sup>*<sup>2</sup> *<sup>T</sup>*<sup>2</sup> and *<sup>η</sup>*<sup>2</sup> *<sup>T</sup>*3, respectively. On the other hand, the terms *Ci*,*<sup>j</sup> Tl*(*k*) <sup>∀</sup>*i*, *<sup>j</sup>*, *<sup>l</sup>* <sup>∈</sup> [1, 2, 3], are the regression coefficients for textures *<sup>T</sup>*1(*t*), *<sup>T</sup>*2(*t*) and *<sup>T</sup>*3(*t*), respectively.

Now let us analyze the variances/covariances of the residual terms *η*<sup>2</sup> *Ti* by means of the following Σ matrix form Equation (7):

$$
\Sigma = \begin{pmatrix}
\Sigma\_{T1}^2 & \Upsilon\_{1,2} & \Upsilon\_{1,3} \\
\Upsilon\_{2,1} & \Sigma\_{T2}^2 & \Upsilon\_{2,3} \\
\Upsilon\_{3,1} & \Upsilon\_{3,2} & \Sigma\_{T3}^2
\end{pmatrix} \tag{7}
$$

where Υ1,2 is the covariance between *η*<sup>2</sup> *<sup>T</sup>*<sup>1</sup> and *<sup>η</sup>*<sup>2</sup> *<sup>T</sup>*<sup>2</sup> (defined as <sup>Υ</sup>1,2 <sup>=</sup> *cov*(*η*<sup>2</sup> *<sup>T</sup>*1, *<sup>η</sup>*<sup>2</sup> *<sup>T</sup>*2)); <sup>Υ</sup>1,3 is the covariance between *η*<sup>2</sup> *<sup>T</sup>*<sup>1</sup> and *<sup>η</sup>*<sup>2</sup> *<sup>T</sup>*<sup>3</sup> (defined as <sup>Υ</sup>1,3 <sup>=</sup> *cov*(*η*<sup>2</sup> *<sup>T</sup>*1, *<sup>η</sup>*<sup>2</sup> *<sup>T</sup>*3)), and so on.

Based on the earlier conditions and using the concept of statistical independence between two random processes at the same time (in pairs), causality can be defined in time. An example of the causality between *T*1 and *T*2 is as in the following expression:

$$F\_{T2,T1} = \ln\left[\frac{\Sigma\_{T1}^1 \times \Sigma\_{T2}^1}{\Sigma\_{T1^2 \times \Sigma\_{T2}^2}}\right] \tag{8}$$

The Equation (8) is commonly known as the causality in the time domain. From this equation, if the random processes *<sup>T</sup>*1(*t*) and *<sup>T</sup>*2(*t*) are statistically independent, then *FT*1,*T*<sup>2</sup> = 0; otherwise there will be causality from one to another.

In the Equation (1), Σ<sup>1</sup> *<sup>T</sup>*<sup>1</sup> measures the precision of the autoregressive model to predict *<sup>T</sup>*1(*t*), established on the past samples.

Then again, Σ<sup>2</sup> *<sup>T</sup>*<sup>1</sup> in the expression (4) measures the precision to predict *<sup>T</sup>*1(*t*) based on the previous values of *T*1(*t*), *T*2(*t*) and *T*3(*t*) at the same time. Returning to the case of taking only 2 textures at the same time *T*1(*t*) and *T*2(*t*) and according to References [37] and [7], if Σ<sup>2</sup> *<sup>T</sup>*<sup>2</sup> < <sup>Σ</sup><sup>1</sup> *<sup>T</sup>*<sup>1</sup> then it is said that *<sup>T</sup>*2(*t*) has a causal influence on *T*1(*t*). The causality is defined by the following equation:

$$F\_{T2\rightarrow T1} = \ln\left[\frac{\Sigma\_{T1}^1}{\Sigma\_{T1}^2}\right] \tag{9}$$

It is relatively easy to see that if *FT*<sup>2</sup>→*T*<sup>1</sup> <sup>=</sup> 0 then there is no causal influence from *<sup>T</sup>*2(*t*) towards *T*1(*t*), at any other values, the result will be otherwise. On the other hand, the causal influence of *T*1(*t*) towards *T*2(*t*) is established using the following equation:

$$F\_{T1 \to T2} = \ln \left[ \frac{\Sigma\_{T2}^1}{\Sigma\_{T2}^2} \right] \tag{10}$$

### **3. Methodology**

In the current section, we describe the methodology developed for the WGC technique with a GA support applied to natural scenery.

For the use of the CBIR, there are different determining factors that must be taken into account while extracting the information from the images, such as luminosity, orientation, scale, homogeneity, and so forth. The main characteristic in our proposed patterns is texture, such that we try to create a base dictionary to later create the time series from the reading of the images and their comparison with the dictionary, with which the theory of WGC was applied.

For the development of the dictionary, a set of *k* textures are manually selected on the images to be studied, which we will call *reference textures*. The *k* generated textures represent parts of objects such as the sky, clouds, grass, rock, and so forth, trying to make a manual segmentation of the scenery as shown in Figure 3. In Section 6, the *k* = 4 textures test is shown for 6 scenery-classes.

Once the set of the *k* reference textures has been obtained, the values in the HSI color space of each of them are examined to create a range of maximum and minimum values which represent them, these values help us to define the thresholds of comparison for the test textures of a query image.

**Figure 3.** Example of segmentation texture zone in a natural scene.

The proposed methodology for the identification and classification of scenery by WGC is shown in Figure 4. The blocks of the architecture are described below.

	- Grid image. The work done in Reference [21] is taken as a reference, a regular grid of 10 × 10 windows is considered for the CBIR texture analysis; in our proposal, we use a grid of *r* × *c* windows, which has the property of *<sup>r</sup>* <sup>=</sup> *<sup>c</sup>*, where *<sup>c</sup>* :=*number of windows in horizontal (columns)*, and *r* :=*number of windows in vertical (lines)*.


**Figure 4.** Learning and testing architecture for the classification system.

6. Generation of time series.For each *M<sup>i</sup> <sup>s</sup>* of the previous step, each matrix entry is compared to the *<sup>k</sup>*−textures of the dictionary to construct a discrete signal as a time series *<sup>T</sup><sup>i</sup> <sup>S</sup>*, defined as a matrix of size *k* × *w*.

In comparison, the value 1 is assigned if the feature neighborhood approaches the dictionary texture and 0 if not, according to the threshold values which characterize each texture as they were previously presented. After processing the entries of all *M<sup>i</sup> <sup>s</sup>*, the set of signals for each scenery is stored in the *FMs*, the time series matrix corresponding to a class *s* that contains *Imgs* images.

7. Wiener-Granger Causality analysis. Each *FMs* matrix created in the previous step was carried to the WGC analysis to obtain the causal relationships, contained among each one of the base textures. A matrix of causality relationships, η*s*, related to the training images was generated, as shown in Figure 7; therein, darker colors represent stronger relationships and these can be depicted through a state diagram where continuous lines represent only the stronger ones. The analysis of causality was computed with the causality toolbox MVGC [38], which was invoked as an external system call.

Once the causality analysis has been made for each of the *Cs* scenery, we get a causality relationships matrix η*<sup>s</sup>* of size *k* × *k*, with the total of the causal relationships *FTi*,*Tj* from the texture *Ti* → *Tj* (as given in Equation (11)), such that if a value of *FTi*,*Tj* <sup>=</sup> 0 means that there is no causal relationship of the texture *i* → *j*, and in the measure that the value increases with respect to other η*<sup>s</sup>* values, we say that the causal relationship is significant with respect to others.

$$
\eta\_{\mathfrak{k}} = \begin{bmatrix} F\_{T\_1, T\_1} & F\_{T\_1, T\_2} & \dots & F\_{T\_1, T\_k} \\ F\_{T\_2, T\_1} & F\_{T\_2, T\_2} & \dots & F\_{T\_2, T\_k} \\ \vdots & \vdots & \vdots & \vdots \\ F\_{T\_k, T\_1} & F\_{T\_k, T\_2} & \dots & F\_{T\_k, T\_k} \end{bmatrix} \tag{11}
$$

The causality matrices η*<sup>s</sup>* are normalized according to the total sum of their values, being *Ns* = ∑*<sup>k</sup> <sup>i</sup>*,*j*=<sup>1</sup> *FTi*,*Tj*, such that η*<sup>N</sup> <sup>s</sup>* is the normalized matrix of the *<sup>s</sup>*−th scenario, for *<sup>s</sup>* <sup>=</sup> 1, ... , *Cs*, with *Cs*: the number of scenery types considered, as given in the Equation (12). From this resulting matrix the values of the main diagonal are not taken into account because these values do not generate force in the causality relationship; as observed in the theory, there is no causal relationship between the same variables.

At the end, for *Cs* classes or scenery, the *total concentration of the matrices*, Γ, is defined as given in the Equation (13).

*FTk*,*T*<sup>1</sup> *FTk*,*T*<sup>2</sup> ... *FTk*,*Tk*

$$
\mathfrak{m}\_s^N = \begin{bmatrix} F\_{T\_1, T\_1} & F\_{T\_1, T\_2} & \dots & F\_{T\_1, T\_k} \\ F\_{T\_2, T\_1} & F\_{T\_2, T\_2} & \dots & F\_{T\_2, T\_k} \\ \vdots & \vdots & \vdots & \vdots \\ \dots & \dots & \dots & \dots \end{bmatrix} \ast \frac{1}{N\_s} \tag{12}
$$

$$\Gamma = \cup\_{l=1}^{\mathbb{C}s} \mathfrak{n}\_l^N = \{ \mathfrak{n}\_1^N, \mathfrak{n}\_2^N, \dots \mathfrak{n}\_{\mathbb{C}s}^N \}\tag{13}$$

The Γ matrices as entries serve as a descriptive pattern for each scenery or class contained in the database. 8. Selection of causal relationships by means of Genetic Algorithm. To look for the causal relationships among different variables that are more important or relevant, for each of the scenes, this can be accomplished in a simple way by eliminating the relationships that have a numerical value less than a previously established threshold.

However, one disadvantage of this method is the establishment of the threshold to be used, because there is no *a priori* knowledge of the optimal value; in addition, the complexity increases when the number of textures increases in the dictionary, along with the number of classes and images to be examined. Other drawback of this solution is that some of the weak relationships could also be important in order to characterize a scenery. So there is a need to implement an automatic selection which discriminates the relevant relationships as a combinatorial optimization process. Genetic algorithms (GA) have been used successfully in several computer vision problems together with the digital image processing [39] and classification [1,9,40–42]. In this work the GA is also the right solution for the required optimization.

**Figure 5.** A 10 × 10 grid partition image example, every grid has a window of 10 × 10 pixels size.

**Figure 6.** Reading image example among the grid neighborhood of the image.

**Figure 7.** Generation of a texture-based causality relationship matrix, η*s*, using the WGC analysis.

### **4. Genetic Algorithm Proposal**

Looking for the analysis of the Γ matrices generated by the WGC to find the significant causality relationships for one scenery, we propose each matrix to be treated with a GA implementation. In this section we provide the GA proposal in detail.

In this approach, each matrix η*<sup>N</sup> <sup>s</sup>* ∈ Γ is expressed using vector representation, see Figure 8 parts *(a)* and *(b)*; this is achieved only by concatenating the rows of the matrix η*<sup>N</sup> <sup>s</sup>* , then the entries of the diagonals are eliminated as shown in Figure 8 part *(c)*. In Figure 8, part *(d)*, a reallocation of the values after the previous elimination is adjusted. This provides a vector of continuous index having the size of each vector 1 <sup>×</sup> (*k*<sup>2</sup> <sup>−</sup> *<sup>k</sup>*) for each *<sup>s</sup>*−th row, one per scenery. Following this process, finally, the matrix <sup>τ</sup> is

created, which contains the linear conformation of each matrix η*<sup>N</sup> <sup>s</sup>* , with *<sup>s</sup>* <sup>=</sup> {1, 2, ... , *Cs*} in different rows, as shown in Figure 8 part *(e)*.

**Figure 8.** The τ matrix generation process, for every *s*−th row ∈ τ.

### *4.1. Individual Codification*

An individual binary representation for one scenery, *τ*[*i*], is conformed in order to create a filter type array of size 1 <sup>×</sup> (*k*<sup>2</sup> <sup>−</sup> *<sup>k</sup>*) of zeros and ones, such that if an input or causal relationship is selected in that array, the value 1 is used, and 0 if not. So we have *Cs* rows, one row per scenery, it is intended that each row of the filter matrix could be different from the other lines, with the purpose of characterizing each type of scenery in a unique way.

It is then necessary to apply an automatic process to determine which values of the matrix τ are relevant features to distinguish the causal relationships of each scenery, and based on this result, it selects which values are going to be removed for the preset number of textures. With the selection of the most relevant causal values, it is sought to have a classification by means of a distance classifier, towards the matrix τ for each one of the query images.

### *4.2. Fitness Function*

The fitness evaluation of each individual is generated in several parts. First, the Equation (14) is applied to the individual *Gx*, representing a texture relationships selection for the *s*−scenery in question, using the matrix τ in Figure 9.

*<sup>ρ</sup>Gx <sup>s</sup>* <sup>=</sup> *<sup>k</sup>*2−*<sup>k</sup>* ∏ *l*=1 *Gx*(*l*) <sup>∗</sup> <sup>τ</sup>*s*,*<sup>l</sup> Cs* ∑ *m*=1 *<sup>k</sup>*2−*<sup>k</sup>* ∏ *l*=1 *Gx*(*l*) <sup>∗</sup> <sup>τ</sup>*m*,*<sup>l</sup>* , such that *Gx*(*l*) <sup>=</sup> 0 (14)

**Figure 9.** The *Gx* genome construction.

There, <sup>∏</sup> *Gx*(*l*) <sup>∗</sup> <sup>τ</sup>*s*,*<sup>l</sup>* refers to the product of the <sup>τ</sup> entries located at *<sup>s</sup>*−scenery (row *<sup>s</sup>*) and column *<sup>l</sup>*, specifying a causal relationship, accomplishing *Gx*(*l*) is a valid non zero entry of the genome. Thus *<sup>ρ</sup>Gx <sup>s</sup>* is the total probability for the individual *Gx* applied to all scenery.

Based on these data, by means of the probability theory, the individual *Gx* is required to meet the condition: *ρGx <sup>s</sup>* > *ρGx <sup>j</sup>* , such that *<sup>s</sup>* ∈ {1, 2, . . . , *Cs*}, *<sup>s</sup>* <sup>=</sup> *<sup>j</sup>*, and 1 <sup>≤</sup> *<sup>j</sup>* <sup>≤</sup> *Cs*.

That is, *Cs* probabilities corresponding to each scenery evaluating the individual *Gx* are obtained with the calculation of *ρGx <sup>j</sup>* . Equation (15) gives the first step for the optimization process, considering the maximum probability related to the casual relationships which best characterized the *s*−scenery versus the others.

$$\begin{aligned} \bigvee\_{s}^{\mathcal{G}\_{\mathbf{x}}} = \begin{cases} \rho\_{\mathbf{s}}^{\mathcal{G}\_{\mathbf{x}}} & \text{if} \quad \rho\_{\mathbf{s}}^{\mathcal{G}\_{\mathbf{x}}} = \text{Max}\{\rho\_{j}^{\mathcal{G}\_{\mathbf{x}}}\}\_{j=1,\ldots,\mathcal{C}s} \\ 0 & \text{if} \quad \rho\_{\mathbf{s}}^{\mathcal{G}\_{\mathbf{x}}} \neq \text{Max}\{\rho\_{j}^{\mathcal{G}\_{\mathbf{x}}}\}\_{j=1,\ldots,\mathcal{C}s} \end{cases} \end{aligned} \tag{15}$$

Then, the fitness function, *fs*(*Gx*), is determined as Equation (16).

$$f\_{\sf s} \left( \mathbf{G}\_{\sf x} \right) = \begin{cases} \sf C \mathbf{P}\_{\sf s} & \text{if } \begin{array}{c} \sf Gx \\ \sf \\ s \\ 0 \end{array} \\ \sf 0 & \text{if } \begin{array}{c} \sf Gx \\ \sf \\ s \end{array} \\ \sf \end{cases} \tag{16}$$

To this end, the images contained in the *s*−scenery are consulted, using the re-substitution test. Each image query gives the scenery which belongs to filling the information of a confusion matrix that is used to calculate the percentage of classification. The image consult query process is described in the following paragraph. Later, in Section 4.3 the *global fitness* is taken into account for the population evolution in the GA loop process.

#### 4.2.1. Creating a Query from a Single Image

In order to classify an s-scenery image considering the relationships specified in a *Gx* individual, a related causal relationship matrix needs to be constructed.

The first step consists of creating a set of *M* synthetic images, *L*1, *L*2, ... , *LM*, from a single *L* image is performed. This is produced by means of manipulating the first reading of the image, making a circular shift of *d* positions for each new synthetic image, in order to create several samples of the same image as shown in Figure 10. In this way, the respective query matrix of size <sup>|</sup>*<sup>k</sup>* <sup>×</sup> *<sup>w</sup>* <sup>×</sup> *<sup>M</sup>*<sup>|</sup> (*<sup>k</sup>* :<sup>=</sup> number of textures in the dictionary, *w* := number of neighborhoods, and *M* := number of synthetic images) is generated to

feed the WGC analysis process and to obtain the resulting normalized causal relationship matrix η*<sup>N</sup> <sup>L</sup>* of size |*k* × *k*|. These steps are carried out by Equations (11) and (12).

**Figure 10.** Representation of a query image construction into *N* samples.

Then the manipulation of η*<sup>N</sup> <sup>L</sup>* is performed as shown in the stage presented in Figure 8 to obtain the linear representation of the matrix. The last query step consists of applying the *k*-NN classifier (with *k* := 1) to determine which τ scenery (line) has the closest relationship to the linear relationship representation of image *L*, considering only the relationship indicated with the *Gx* non-zero values.

### *4.3. GA Implementation*

A genetic algorithm is applied for each τ line to automatically select the most representative causal relationship of each scenery. Figure 11a shows the general algorithm flowchart of this approach.

An initial population, *PG*, of *sizeP* individuals is randomly generated, where *sizeP* is an odd number and each individual is of size *<sup>k</sup>*<sup>2</sup> − *<sup>k</sup>*, the size in columns of the matrix <sup>τ</sup>.

Then the *PG* individuals are evaluated with the fitness function, Equation (16), for a particular *s*−scenery, considering the total set of images that conform it, as shown in Figure 11b. The individual's fitness is stored inside a fitness array {•}, as in Equation (17). The {•} array is consequently ordered, from highest to lowest, to find the best individual with the highest fitness.

$$\{\bullet\} = \{f\_1(\mathcal{G}\_\mathbf{x}), f\_2(\mathcal{G}\_\mathbf{x}), \dots, f\_{\text{sizeP}}(\mathcal{G}\_\mathbf{x})\},\tag{17}$$

such that *fp*(*Gx*) <sup>≥</sup> 0 for 1 <sup>≤</sup> *<sup>p</sup>* <sup>≤</sup> *sizeP*. For this proposal, size population *sizeP* <sup>=</sup> 21, the genome length is 12, and the number of iterations was *maxGen* = 100 generations.

To generate the new population, (*sizeP* <sup>−</sup> <sup>1</sup>) \$ 2 triplets of random numbers are generated, e.g., {1,5,1} or {2,4,0}, where the first two numbers are the selected individual numbers that generate the new individuals, the third element of the triplet is one of the two possible operations to be executed; either crossover "1" or mutation "0".

*Electronics* **2019**, *8*, 726

**Figure 11.** Flowchart showing the general GA algorithm implementation. (**a**) Implementation of the proposed General flowchart GA, and (**b**) GA implementation for a particular s-scenery.

For derivatives of each triplet, we generate two new individuals if it is by crossover operator, or if it is by mutation operator the two selected individuals are altered separately, in order to have two new elements for the new generation, 1 mutated individual from one single individual.

The crossover operator is applied at one uniform random point of the two participating chromosomes, and the mutation operator is performed over 10% of the elements of a chromosome, as shown in Figure 12.

**Figure 12.** Genetic operators application.

The genetic operations of mutation and crossover are applied to 30% and 70% of the population respectively, favoring the selection of the highest fitness individuals to be reproduced. The individual with the best fitness passes to the next generation applying elitism. In this way, the population will evolve towards a selection of relevant causal relationships to be able to characterize each scenery.

The end of the GA or stop criteria is given when reaching the 100% classification percentage or a number of generations is attained.

After the GA is applied *Cs* times, the individuals that contain the most relevant relationships for each scenery are found. Then, the τ matrix is updated and its entries are replaced by a zero value whenever the corresponding individual entries have a zero and they keep their value in other case.

### **5. Parallel Approach**

In this section, a parallel algorithm to speed-up the performance of the proposed WGC methodology is presented. The parallel approach works on a distributed memory architecture using MPI library; that is, there is a set of processes without shared memory, and these processes work in parallel, and the communication goes through message exchange to determine the relevant causal relationships of all the scenery. Each process can access the NSDB to extract and work with the corresponding set of images. The algorithm complexity in this proposal is given for the Equation (18).

$$\text{MO}(\text{N}\_{\text{class}} \times k \times \text{CIP} \times r \times c \times n \,\text{Im}\,\text{g} \times \text{t}\_{\text{comp}} \times \text{WGC}\_{tb}) \tag{18}$$

where, *Nclass* := number of classes, *<sup>k</sup>* := number of textures, *CIP* := constant for every image processing, *r* := number of rows in the grid, *c* := number of cols in the grid, *nImg* := total number of images, *tcomp* := comparison time against the base textures, *WGCtb* := causality analysis time. (e.g., for an image of size <sup>640</sup>*x*480 pixels, *<sup>r</sup>* = (640/20), *<sup>c</sup>* = (480/30), and *<sup>p</sup>* is the size of the neighborhood, such that 10 <sup>×</sup> <sup>10</sup> implies *p* = 10) That means, if the number of rows *r*, cols *c* in the grid increases, then the number of images *nImag* in NSDB increases, and the number of base textures *k* in the dictionary increases, and the number of classes *Nclass* increases, thus the computational cost increases. In this way it is necessary to conceive a parallel architecture to solve this problem in a large number of images related to Big Data problems.

The Algorithm 1 shows the procedure which is executed simultaneously by each process and the general process of this parallel proposal is depicted in Figure 13.

**Figure 13.** The proposed parallel algorithm structure.

At the beginning (line 2 of Algorithm 1, Figure 13, tag (1)), each process determines the amount set of images to be read (*ImgBlock*), taking into consideration the total number of images (*Total*\_*IMGs*), the total number of processes (*Total*\_*procs*) and the process identifier (*rank*). A single process can work with images belonging to different scenery (e.g., *Total*\_*IMGs* = 700, *Total*\_*procs* = 70, *imgBlock* = 700/10 = 10 for the *NSDB*).

**Algorithm 1** Parallel algorithm for the causality matrix construction.

1: **procedure** CAUSALITY MATRIX CONSTRUCTION(rank) 2: initialization(ImgBlock,Total\_IMGs,Total\_procs,rank); 3: **for** every *i* in ImgBlock **do** 4: *Img*\_*RGBread*(*image*,*s*, *i*); 5: *Img*\_*preprocessing*(*image*); 6: *RGB*\_*to*\_*HSI*(*image*); 7: *M<sup>i</sup> <sup>s</sup>*= Feature\_extraction(image); 8: *F<sup>i</sup> s*= time\_series\_construction(*M<sup>i</sup> <sup>s</sup>*, *Texture*\_*Dictionary*); 9: Save\_TimeSeries(*F<sup>i</sup> s*); 10: **end for** 11: *Barrier*\_*synchronization*(); 12: **if** rank in {Scenery coordinator ranks} **then** 13: *Fs*= Load\_all\_time\_series(*s*); 14: η*s*= *System*\_*call*(*MVGC*(*Fs*)); 15: τ*s*= Fitting(η*s*); 16: Genetic\_Algorithm(τ*s*); 17: Send(τ*s*,*s*, *General*\_*coordinator rank*); 18: **end if** 19: **if** (rank == *General*\_*coordinator rank*) **then** 20: **for** (every *id* in {Scenery coordinator ranks}) **do** 21: Recv(τ*s*,*s*, *id*); 22: τ(*s*)= τ*s*; 23: **end for** 24: **end if** 25: **end procedure**

Each process works simultaneously with the section of the NSDB, *ImgBlock*, which was assigned to it, performing the following steps (lines 3–10). The reading of the *i*-image in RGB space is the first action to be executed, next up the scenery, *s*, such that *i*-image ∈ *s*, is also obtained (line 4). Then the image preprocessing and the conversion from RGB to HSI domains are carried out (lines 5 and 6, respectively). In line 7, the statistical features are calculated, including the construction of the image grid and neighborhoods, then the CBIR features per each neighborhood generates the *M<sup>i</sup> <sup>s</sup>* matrix; Figure 13, tag (2), represents the execution of lines 4 to 7 of Algorithm 1. Then, *M<sup>i</sup> <sup>s</sup>* and the texture dictionary are used to construct the respective time series, *F<sup>i</sup> <sup>s</sup>*, that is stored on file (lines 8,9 of the algorithm, Figure 13, tag (3)).

Up to this point, all processes work independently; however, in order to ensure that every process has fully accomplished its task, a parallel barrier synchronization (line 11 of the algorithm, Figure 13, tag (4)) should be introduced before continuing with the next step. Here (line 12), only the processes identified as *scenery coordinators* (one process per scenery) continue with the construction of the corresponding *Fs* matrix (line 13, Figure 13, tag (5)), by loading the respective set of *F<sup>i</sup> <sup>s</sup>* matrices (one per scenery image), previously generated. Then, a system call is performed (line 14, Figure 13, tag (6)) to run the MVGC toolbox and obtain the causality relationship matrix, η*s*, from the WGC analysis.

The *Fit*(η*s*) function in line 15 (Figure 13, tag (7)) is in charge of normalization and vector representation of the causality relationship matrix, η*s*. The respective τ*<sup>s</sup>* is thus generated, corresponding to the *s*-th row (scenery) of τ matrix. Line 16 of Algorithm 1 (Figure 13, tag (8)) shows the GA call that is executed by each one of the *Cs* scenery coordinator processes, with *Cs* being the number of scenery. After identifying the most relevant causal relationships by means of the GA, τ*<sup>s</sup>* is updated and sent to the general coordinator (line 17) through a message.

Finally, in lines 19–24, the general coordinator process receives, by means of several messages, the results generated by the scenery coordinator processes (Figure 13, tag (9)). When all message receptions are achieved, the matrix τ is successfully constructed.

### **6. Experimental Results**

The proposal evaluation was generated using the computer power of a 19-processor dual core cluster. Each processor is an Intelc Xeonc CPU E5-2670 v3 2.30 GHz, and 74 GB RAM.

Four image textures *k* = 4 where selected to conform the base dictionary, as shown in Table 1. For each texture the generated values were obtained manually within the images of the database, a set of 20 texture samples were taken from a set of 5 images per class, from each texture the average was extracted in the layer H plus twice the standard deviation, with this the maximum and minimum threshold values for each texture were generated.


**Table 1.** HSI ranges of the base texture dictionary.

The *NSDB* used for the evaluation consists of the following data:


The images were adapted so that some typical classification challenges were considered. The whole set of images was tested in a normal state and introducing Gaussian noise (GN), salt and pepper noise (S&P) of 1%, 3%, and 10% levels respectively, as shown in Figure 14. A rotation transformation was also introduced on each image considering 0◦, 45◦, 90◦, 135◦, and 180◦, as shown in Figure 15. An image consult query was performed following the same procedure described in Section 4.2.1.

**Figure 14.** Example of images with Gaussian, salt, and pepper noise of 1%, 3%, and 10%, respectively.

The results in this section are organized as follows. The image classification performance obtained when applying the WGC theory is first presented. Then the execution times of the proposed parallel methodology are shown.

**Figure 15.** Rotation degrees applied to the images set.

#### *6.1. Classification Results*

To show that the proposed GA implementation was a good solution to select some relevant texture relationships describing a scenery, we compared our proposal (GA version) to the manual strategy (Manual version) introduced in Reference [36]; under the manual strategy only the highest relationship values were selected, establishing a specific threshold. In both versions, the methodology presented in Section 3 for the construction of τ matrix, was executed. Table 2 shows the resulting τ values.

Because there were no *a priori* criteria to determine a threshold value, in the Manual version 25% of the less significant causal relationships per scenery were deleted. Table 3 shows the updated τ matrix after the manual selection.


**Table 2.** The obtained τ matrix values.

**Table 3.** The τ matrix resulting from the manual selection of the highest significant causal relationships.


When executing the GA version for the selection of τ matrix relationships, a larger space of solutions was explored trying to look for the causal relationships that best represent one scenery. The obtained individuals are presented in Table 4, and the updated τ matrix is shown in Table 5.

Both, GA and manual versions where tested using 300 images, 50 per scenery. The manual version only obtained obtaining a 12.53% general classification percentage. The confusion matrix showing the image association per scenery can be seen in Table 6; we observe that most of the images were associated to the coast scenery, and as a result the manual selection test gave a poor classification percentage.


**Table 4.** The best individuals resulting from the evaluation of the GA per scenery.

**Table 5.** Final τ values, applying the better individuals of the GA.


**Table 6.** Confusion matrix for a test with 50 images per scenery, using a manual selection of causal relationships.


With the information in the Table 5, the most representative relations of each natural scenery are generated as a visual representation, the graphs representing the intensity of the causal relations between the *k* = 4 base textures of the dictionary. These graphs will show how textures are related within the corresponding scenery, obtaining the pattern which represents each of them, as it can be appreciated in Figure 16.

**Figure 16.** Evolutionary texture causal relationship graphs resulted for each scenery.

Given these first results, it can be observed that not necessarily the relationships with higher values were the best ones to be selected.

To measure the technical efficiency of our proposal using the GA, the *Recall* (managed as classification percentage), *Precision*, *Accuracy* and *F1 Score* were estimated from the confusion matrices of every test.

The classification results of Figure 17 show that that rotating the images by 45◦, 90◦ and 315◦, the classification performance decays significantly. Also, the noise (GN and S&P) significantly alters the classification which is expected in natural scenery images since the texture is a representative of the type of

image, and the alterations with noise on it degenerate into another possible meaning. In normal conditions, avoiding noise and rotations, the classification performance reaches 100%.

Additionally, Figure 18 shows the estimations of the precision (Figure 18a), recall (Figure 18b), accuracy (Figure 18c), and F1 Score (Figure 18d) averages for the classes contained in the NSDB. In general, the classification of ideal images (0◦) without rotations and noise obtains 100% classification, However, when rotations and noise are added this percentage decreases, particularly the sky class is the most affected.

**Figure 17.** Classification results of our proposal using the GA, considering different noise type and rotation configurations.

**Figure 18.** Technical efficiency measures for the best GA individuals for each scenery: (**a**) precision measure, (**b**) recall measure, (**c**) accuracy measure and (**d**) F1 score measure.

The average result for the GA evaluation is depicted in Figure 19. Figure 19a shows the fitness evolution within a run with 100 generations. Figure 19b shows that, while in some classes the highest fitness is achieved in the first iteration, in the other ones 100 generations are not enough to achieve the best fitness. Figure 19c shows the best fitness obtained through 200 runs considering 100 generations per run and population size set to 21 individuals; all fitness converges near the expected value of 100% classification.

(**a**)Representation of a single execution of the GA. (**b**) Convergence of a simple execution of the GA.

(**c**) Evaluation of the GA with 200 executions and 100 generations per execution.

**Figure 19.** Performance of the GA for the natural scenery contained in the NSDB.

### *6.2. Parallel Methodology Performance*

The execution time taken by the proposed parallel causality methodology applied to the identification and classification of natural scenery, for a total of 700 images, 6 different scenery, and varying the number of processes, are shown in Figure 20. These values were the average time taken for 200 executions. We can observe that the execution time decreased rapidly while increasing the number of processes, getting the best execution time when 125 concurrent processes were defined. With this configuration each process worked with 5 or 6 images, favoring the internal scheduling that used more efficiently the computer resources so that the sequential version execution time decreased by 88.9%.

**Figure 20.** Performance of the parallel methodology algorithm.

### **7. Discussion**

To find the texture causality relationships for characterizing a natural scenery, we found the necessity of GA implementation. With this solution, the automatic discrimination process for selecting the causal relationships that are important or relevant for the classification of the proposed scenery was successfully achieved. Compared to some of the articles in the literature, the possibility of our approach to perform the selection in an automatic way allows the study of the scenery classification problem considering more parameters. In this way, a larger number of scenery or base textures for future implementations of several purposes of texture classification could be taken into account considering the efficient methodology and evolutionary algorithm proposed in this work. This proposal could help in recognizing natural scenarios in the navigation of an autonomous car or possibly a drone, being an important element in the safety of autonomous vehicles navigation.

As we can see in Figure 17, the classification percentage obtained by means of the selected features for the evolutionary process surpasses those obtained by the manual selection version. This result is important because the representative causal relationships of the scenery are selected in such a way that they numerically escape the manual perspective; that is to say, in the manual selection strategy a non-significant threshold value is specified, such that any value lower than the threshold is set to zero, but evolutionary strategy turns out that some of these causal relationships are relevant to classify the scenery, marking differences with another similar scenery. From the causality theory applied to the image reading sequences, we are trying to infer the order of appearance of the textures typified in the base dictionary seeking to represent them as temporal visual reading that we see as a type of natural scenery.

### **8. Conclusions**

In this paper a novel proposal for the use of the Wiener-Granger causality theory supported by a genetic algorithm was presented, along with the CBIR self-content analysis, applied for the identification and classification of 6 natural scenery: coast, forest, mountain, prairie, river/lake, and sky/cloud. Considering the new formulation it was possible to find a set of descriptors from the causality matrices in order to represent a scenery class, from a base set of reference textures, proposing a characterization of images based on the continuous appearance of textures within them; the base dictionary in this approach included the textures: Cloud, Sky, Rock, and Forest. Unlike others approaches, our methodology deals with the rotation and image noise considerations, and the results show excellent classification percentages.

Under this approach we have 100% image classification for the whole dataset, and the methodology provided the next good classification rate for 180◦ rotation, and the sensitivity for intermediate rotation levels (45◦, 90◦), and had good results for the salt and pepper image noise.

In relation to the proposal performance, the design of a parallel computing algorithm was developed. A reduction in execution times was achieved using a 19-processor dual-core cluster server, and the MPI tool, reaching an 88.9% decrease of the sequential version execution time when 125 processes were launched.

Future work includes the study of performance of this proposal using other parallel architectures; e.g., the GPU technology could perform efficiently for the image feature extraction stage, as well as the implementation of other evolutionary algorithms, such as Genetic Programming in order to analyze all together the image textures looking to characterize the whole scenery and its associations with paradigm of visual comprehension.

**Author Contributions:** Writing–review and editing, C.B.-A.; investigation, C.B.-A., C.A.-C. and J.V.-C.; resources, J.V.-C. and C.A.-C.; writing–original draft preparation, C.B.-A; validation, G.R.-A., C.A.-C. and J.V.-C.; conceptualization G.R.-A, C.B.-A., C.A.-C.; formal analysis, C.A.-C., G.R.-A., J.V.-C.; methodology, C.B.-A., J.V.-C. and G.R.-A; C.B.-A. supervised the overall research work. All authors contributed to the discussion and conclusion of this research.

**Funding:** This research received no external funding.

**Acknowledgments:** Cesar Benavides thanks the CONACyT for the scholarship support. This work has been supported by Fundación Carolina, Spain, under the scholarship program 2016–2017. This work has been done under project EL006-18, granted by the Metropolitan Autonomous University, Unidad Azcapotzalco, Mexico.

**Conflicts of Interest:** The authors declare that there is no conflict of interests regarding the publication of this paper.

### **References**




c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
