Multithreshold Segmentation by Using an Algorithm Based on the Behavior of Locust Swarms

Abstract

As an alternative to classical techniques, the problem of image segmentation has also been handled through evolutionary methods. Recently, several algorithms based on evolutionary principles have been successfully applied to image segmentation with interesting performances. However, most of them maintain two important limitations: (1) they frequently obtain suboptimal results (misclassifications) as a consequence of an inappropriate balance between exploration and exploitation in their search strategies; (2) the number of classes is fixed and known in advance. This paper presents an algorithm for the automatic selection of pixel classes for image segmentation. The proposed method combines a novel evolutionary method with the definition of a new objective function that appropriately evaluates the segmentation quality with respect to the number of classes. The new evolutionary algorithm, called Locust Search (LS), is based on the behavior of swarms of locusts. Different to the most of existent evolutionary algorithms, it explicitly avoids the concentration of individuals in the best positions, avoiding critical flaws such as the premature convergence to suboptimal solutions and the limited exploration-exploitation balance. Experimental tests over several benchmark functions and images validate the efficiency of the proposed technique with regard to accuracy and robustness.

Erik Cuevas 1, 2 and Adrian Gonzalez 1 and Fernando Fausto 1 and Daniel Zaldivar 1, 2 and Marco Perez-Cisneros 3

Academic Editor:Bogdan Dumitrescu

1, Departamento de Electronica, Universidad de Guadalajara, CUCEI, Avenida Revolucion 1500, 44430 Guadalajara, JAL, Mexico
2, Centro Tapatio Educativo AC, Vicente Guerrero 232, Colonia Zapopan Centro, 45100 Zapopan, JAL, Mexico
3, Centro Universitario de Tonala, Avenida Nuevo Periferico No. 555 Ejido San Jose Tatepozco, 48525 Tonala, JAL, Mexico

Received 7 October 2014; Revised 25 May 2015; Accepted 10 June 2015; 27 August 2015

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Image segmentation [1] consists in grouping image pixels based on some criteria such as intensity, color, and texture and still represents a challenging problem within the field of image processing. Edge detection [2], region-based segmentation [3], and thresholding methods [4] are the most popular solutions for image segmentation problems.

Among such algorithms, thresholding is the simplest method. It works by considering threshold (points) values to adequately separate distinct pixels regions within the image being processed. In general, thresholding methods are divided into two types depending on the number of threshold values, namely, bilevel and multilevel. In bilevel thresholding, only a threshold value is required to separate the two objects of an image (e.g., foreground and background). On the other hand, multilevel thresholding divides pixels into more than two homogeneous classes that require several threshold values.

The thresholding methods use a parametric or nonparametric approach [5]. In parametric approaches [6, 7], it is necessary to estimate the parameters of a probability density function that is capable of modelling each class. A nonparametric technique [8-11] employs a given criteria such as the between-class variance or the entropy and error rate, in order to determine optimal threshold values.

A common method to accomplish parametric thresholding is the modeling of the image histogram through a Gaussian mixture model [12] whose parameters define a set of pixel classes (threshold points). Therefore, each pixel that belongs to a determined class is labeled according to its corresponding threshold points with several pixel groups gathering those pixels that share a homogeneous grayscale level.

The problem of estimating the parameters of a Gaussian mixture that better model an image histogram has been commonly solved through the Expectation Maximization (EM) algorithm [13, 14] or gradient-based methods such as Levenberg-Marquardt, LM [15]. Unfortunately, EM algorithms are very sensitive to the choice of the initial values [16], meanwhile gradient-based methods are computationally expensive and may easily get stuck within local minima [17].

As an alternative to classical techniques, the problem of Gaussian mixture identification has also been handled through evolutionary methods. In general, they have demonstrated to deliver better results than those based on classical approaches in terms of accuracy and robustness [18]. Under these methods, an individual is represented by a candidate Gaussian mixture model. Just as the evolution process unfolds, a set of evolutionary operators are applied in order to produce better individuals. The quality of each candidate solution is evaluated through an objective function whose final result represents the similarity between the mixture model and the histogram. Some examples of these approaches involve optimization methods such as Artificial Bee Colony (ABC) [19], Artificial Immune Systems (AIS) [20], Differential Evolution (DE) [21], Electromagnetism Optimization (EO) [22], Harmony Search (HS) [23], and Learning Automata (LA) [24]. Although these algorithms own interesting results, they present two important limitations. (1) They frequently obtain suboptimal approximations as a consequence of a limited balance between exploration and exploitation in their search strategies. (2) They are based on the assumption that the number of Gaussians (classes) in the mixture is preknown and fixed; otherwise, they cannot work. The cause of the first limitation is associated with their evolutionary operators employed to modify the individual positions. In such algorithms, during their evolution, the position of each agent for the next iteration is updated yielding an attraction towards the position of the best particle seen so far or towards other promising individuals. Therefore, as the algorithm evolves, these behaviors cause that the entire population rapidly concentrates around the best particles, favoring the premature convergence and damaging the appropriate exploration of the search space [25, 26]. The second limitation is produced as a consequence of the objective function that evaluates the similarity between the mixture model and the histogram. Under such an objective function, the number of Gaussians functions in the mixture is fixed. Since the number of threshold values (Gaussian functions) used for image segmentation varies depending on the image, the best threshold number and values are obtained by an exhaustive trial and error procedure.

On the other hand, bioinspired algorithms represent a field of research that is concerned with the use of biology as a metaphor for producing optimization algorithms. Such approaches use our scientific understanding of biological systems as an inspiration that, at some level of abstraction, can be represented as optimization processes.

In the last decade, several optimization algorithms have been developed by a combination of deterministic rules and randomness, mimicking the behavior of natural phenomena. Such methods include the social behavior of bird flocking and fish schooling such as the Particle Swarm Optimization (PSO) algorithm [27] and the emulation of the differential evolution in species such as the Differential Evolution (DE) [28]. Although PSO and DE are the most popular algorithms for solving complex optimization problems, they present serious flaws such as premature convergence and difficulty to overcome local minima [29, 30]. The cause for such problems is associated with the operators that modify individual positions. In such algorithms, during the evolution, the position of each agent for the next iteration is updated yielding an attraction towards the position of the best particle seen so far (in case of PSO) or towards other promising individuals (in case of DE). As the algorithm evolves, these behaviors cause that the entire population rapidly concentrates around the best particles, favoring the premature convergence and damaging the appropriate exploration of the search space [31, 32].

Recently, the collective intelligent behavior of insect or animal groups in nature has attracted the attention of researchers. The intelligent behavior observed in these groups provides survival advantages, where insect aggregations of relatively simple and "unintelligent" individuals can accomplish very complex tasks using only limited local information and simple rules of behavior [33]. Locusts (Schistocerca gregaria ) are a representative example of such collaborative insects [34]. Locust is a kind of grasshopper that can change reversibly between a solitary and a social phase, with clear behavioral differences among both phases [35]. The two phases show many differences regarding the overall level of activity and the degree to which locusts are attracted or repulsed among them [36]. In the solitary phase, locusts avoid contact to each other (locust concentrations). As consequence, they distribute throughout the space, exploring sufficiently over the plantation [36]. On the other hand, in the social phase, locusts frantically concentrate around those elements that have already found good food sources [37]. Under such a behavior, locusts attempt to efficiently find better nutrients by devastating promising areas within the plantation.

This paper presents an algorithm for the automatic selection of pixel classes for image segmentation. The proposed method combines a novel evolutionary method with the definition of a new objective function that appropriately evaluates the segmentation quality with regard to the number of classes. The new evolutionary algorithm, called Locust Search (LS), is based on the behavior presented in swarms of locusts. In the proposed algorithm, individuals emulate a group of locusts which interact to each other based on the biological laws of the cooperative swarm. The algorithm considers two different behaviors: solitary and social. Depending on the behavior, each individual is conducted by a set of evolutionary operators which mimics different cooperative conducts that are typically found in the swarm. Different to most of existent evolutionary algorithms, the behavioral model in the proposed approach explicitly avoids the concentration of individuals in the current best positions. Such fact allows avoiding critical flaws such as the premature convergence to suboptimal solutions and the incorrect exploration-exploitation balance. In order to automatically define the optimal number of pixel classes (Gaussian functions in the mixture), a new objective function has been also incorporated. The new objective function is divided into two parts. The first part evaluates the quality of each candidate solution in terms of its similarity with regard to the image histogram. The second part penalizes the overlapped area among Gaussian functions (classes). Under these circumstances, Gaussian functions that do not "positively" participate in the histogram approximation could be easily eliminated in the final Gaussian mixture model.

In order to illustrate the proficiency and robustness of the proposed approach, several numerical experiments have been conducted. Such experiments are divided into two sections. In the first part, the proposed LS method is compared to other well-known evolutionary techniques over a set of benchmark functions. In the second part, the performance of the proposed segmentation algorithm is compared to other segmentation methods which are also based on evolutionary principles. The results in both cases validate the efficiency of the proposed technique with regard to accuracy and robustness.

This paper is organized as follows: in Section 2 basic biological issues of the algorithm analogy are introduced and explained. Section 3 describes the novel LS algorithm and its characteristics. A numerical study on different benchmark function is presented in Section 4 while Section 5 presents the modelling of an image histogram through a Gaussian mixture. Section 6 exposes the LS segmentation algorithm and Section 7 the performance of the proposed segmentation algorithm. Finally, Section 8 draws some conclusions.

2. Biological Fundamentals and Mathematical Models

Social insect societies are complex cooperative systems that self-organize within a set of constraints. Cooperative groups are good at manipulating and exploiting their environment, defending resources and breeding, yet allowing task specialization among group members [38, 39]. A social insect colony functions as an integrated unit that not only possesses the ability to operate at a distributed manner but also undertakes a huge construction of global projects [40]. It is important to acknowledge that global order for insects can arise as a result of internal interactions among members.

Locusts are a kind of grasshoppers that exhibit two opposite behavioral phases: solitary and social (gregarious). Individuals in the solitary phase avoid contact to each other (locust concentrations). As consequence, they distribute throughout the space while sufficiently exploring the plantation [36]. In contrast, locusts in the gregarious phase gather into several concentrations. Such congregations may contain up to 10¹⁰ members, cover cross-sectional areas of up to 10 km² , and a travelling capacity up to 10 km per day for a period of days or weeks as they feed causing a devastating crop loss [41]. The mechanism to switch from the solitary phase to the gregarious phase is complex and has been a subject of significant biological inquiry. Recently, a set of factors has been implicated to include geometry of the vegetation landscape and the olfactory stimulus [42].

Only few works [36, 37] that mathematically model the locust behavior have been published. Both approaches develop two different minimal models with the goal of reproducing the macroscopic structure and motion of a group of locusts. Considering that the method in [36] focuses on modelling the behavior for each locust in the group, its fundamentals have been employed to develop the algorithm that is proposed in this paper.

2.1. Solitary Phase

This section describes how each locust's position is modified as a result of its behavior under the solitary phase. Considering that [figure omitted; refer to PDF] represents the current position of the [figure omitted; refer to PDF] th locust in a group of [figure omitted; refer to PDF] different elements, the new position [figure omitted; refer to PDF] is calculated by using the following model: [figure omitted; refer to PDF] with [figure omitted; refer to PDF] corresponding to the change of position that is experimented by [figure omitted; refer to PDF] as a consequence of its social interaction with all other elements in the group.

Two locusts in the solitary phase exert forces on each other according to basic biological principles of attraction and repulsion (see, e.g., [36]). Repulsion operates quite strongly over a short length scale in order to avoid concentrations. Attraction is weaker and operates over a longer length scale, providing the social force that is required to maintain the group's cohesion. Therefore, the strength of such social forces can be modeled by the following function: [figure omitted; refer to PDF] Here, [figure omitted; refer to PDF] is a distance, [figure omitted; refer to PDF] describes the strength of attraction, and [figure omitted; refer to PDF] is the typical attractive length scale. We have scaled the time and the space coordinates so that the repulsive strength and length scale are both represented by the unity. We assume that [figure omitted; refer to PDF] and [figure omitted; refer to PDF] so that repulsion is stronger and features in a shorter-scale, while attraction is applied in a weaker and longer-scale; both facts are typical for social organisms [21]. The social force exerted by locust [figure omitted; refer to PDF] over locust [figure omitted; refer to PDF] is [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the distance between the two locusts and [figure omitted; refer to PDF] is the unit vector pointing from [figure omitted; refer to PDF] to [figure omitted; refer to PDF] . The total social force on each locust can be modeled as the superposition of all of the pairwise interactions: [figure omitted; refer to PDF] The change of position [figure omitted; refer to PDF] is modeled as the total social force experimented by [figure omitted; refer to PDF] as the superposition of all of the pairwise interactions. Therefore, [figure omitted; refer to PDF] is defined as follows: [figure omitted; refer to PDF] In order to illustrate the behavioral model under the solitary phase, Figure 1 presents an example, assuming a population of three different members ( [figure omitted; refer to PDF] ) which adopt a determined configuration in the current iteration [figure omitted; refer to PDF] . As a consequence of the social forces, each element suffers an attraction or repulsion to other elements depending on the distance among them. Such forces are represented by [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] . Since [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are too close, the social forces [figure omitted; refer to PDF] and [figure omitted; refer to PDF] present a repulsive nature. On the other hand, as the distances [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are quite long, the social forces [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] between [figure omitted; refer to PDF] and [figure omitted; refer to PDF] all belong to the attractive nature. Therefore, the change of position [figure omitted; refer to PDF] is computed as the vector resultant between [figure omitted; refer to PDF] and [figure omitted; refer to PDF] ( [figure omitted; refer to PDF] ) is [figure omitted; refer to PDF] . The values [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are also calculated accordingly.

Figure 1: Behavioral model under the solitary phase.

[figure omitted; refer to PDF]

In addition to the presented model [36], some studies [43-45] suggest that the social force [figure omitted; refer to PDF] is also affected by the dominance of the involved individuals [figure omitted; refer to PDF] and [figure omitted; refer to PDF] in the pairwise process. Dominance is a property that relatively qualifies the capacity of an individual to survive, in relation to other elements in a group. The locust's dominance is determined by several characteristics such as size, chemical emissions, and location with regard to food sources. Under such circumstances, the social force is magnified or weakened depending on the most dominant individual that is involved in the repulsion-attraction process.

2.2. Social Phase

In this phase, locusts frantically concentrate around the elements that have already found good food sources. They attempt to efficiently find better nutrients by devastating promising areas within the plantation. In order to simulate the social phase, the food quality index [figure omitted; refer to PDF] is assigned to each locust [figure omitted; refer to PDF] of the group as such index reflects the quality of the food source where [figure omitted; refer to PDF] is located.

Under this behavioral model, each of the [figure omitted; refer to PDF] elements of the group is ranked according to its corresponding food quality index. Afterwards, the [figure omitted; refer to PDF] elements featuring the best food quality indexes are selected ( [figure omitted; refer to PDF] ). Considering a concentration radius [figure omitted; refer to PDF] that is created around each selected element, a set of [figure omitted; refer to PDF] new locusts is randomly generated inside [figure omitted; refer to PDF] . As a result, most of the locusts will be concentrated around the best [figure omitted; refer to PDF] elements. Figure 2 shows a simple example of the behavioral model under the social phase. In the example, the configuration includes eight locusts ( [figure omitted; refer to PDF] ), just as it is illustrated by Figure 2(a) that also presents the food quality index for each locust. A food quality index near to one indicates a better food source. Therefore, Figure 2(b) presents the final configuration after the social phase, assuming [figure omitted; refer to PDF] .

Figure 2: Behavioral model under the social phase. (a) Initial configuration and food quality indexes, (b) final configuration after the operation of the social phase.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

3. The Locust Search (LS) Algorithm

In this paper, some behavioral principles drawn from a swarm of locusts have been used as guidelines for developing a new swarm optimization algorithm. The LS assumes that entire search space is a plantation, where all the locusts interact to each other. In the proposed approach, each solution within the search space represents a locust position inside the plantation. Every locust receives a food quality index according to the fitness value of the solution that is symbolized by the locust's position. As it has been previously discussed, the algorithm implements two different behavioral schemes: solitary and social. Depending on each behavioral phase, each individual is conducted by a set of evolutionary operators which mimics the different cooperative operations that are typically found in the swarm. The proposed method formulates the optimization problem in the following form: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is a nonlinear function whereas [figure omitted; refer to PDF] is a bounded feasible search space, which is constrained by the lower ( [figure omitted; refer to PDF] ) and upper ( [figure omitted; refer to PDF] ) limits.

In order to solve the problem formulated in (6), a population [figure omitted; refer to PDF] of [figure omitted; refer to PDF] locusts (individuals) is evolved inside the LS operation from the initial point ( [figure omitted; refer to PDF] ) to a total [figure omitted; refer to PDF] number of iterations ( [figure omitted; refer to PDF] ). Each locust [figure omitted; refer to PDF] ( [figure omitted; refer to PDF] ) represents an [figure omitted; refer to PDF] -dimensional vector [figure omitted; refer to PDF] where each dimension corresponds to a decision variable of the optimization problem to be solved. The set of decision variables constitutes the feasible search space [figure omitted; refer to PDF] , where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] correspond to the lower and upper bounds for the dimension [figure omitted; refer to PDF] , respectively. The food quality index that is associated with each locust [figure omitted; refer to PDF] (candidate solution) is evaluated through an objective function [figure omitted; refer to PDF] whose final result represents the fitness value of [figure omitted; refer to PDF] . In the LS algorithm, each iteration of the evolution process consists of two operators: (A) solitary and (B) social. Beginning by the solitary stage, the set of locusts is operated in order to sufficiently explore the search space. On the other hand, the social operation refines existent solutions within a determined neighborhood (exploitation) subspace.

3.1. Solitary Operation (A)

One of the most interesting features of the proposed method is the use of the solitary operator to modify the current locust positions. Under this approach, locusts are displaced as a consequence of the social forces produced by the positional relations among the elements of the swarm. Therefore, near individuals tend to repel each other, avoiding the concentration of elements in regions. On the other hand, distant individuals tend to attract to each other, maintaining the cohesion of the swarm. A clear difference to the original model in [20] considers that social forces are also magnified or weakened depending on the most dominant (best fitness values) individuals that are involved in the repulsion-attraction process.

In the solitary phase, a new position [figure omitted; refer to PDF] is produced by perturbing the current locust position [figure omitted; refer to PDF] with a change of position [figure omitted; refer to PDF] . The change of position [figure omitted; refer to PDF] is the result of the social interactions experimented by [figure omitted; refer to PDF] as a consequence of its repulsion-attraction behavioral model. Such social interactions are pairwise computed among [figure omitted; refer to PDF] and the other [figure omitted; refer to PDF] individuals in the swarm. In the original model, social forces are calculated by using (3). However, in the proposed method, it is modified to include the best fitness value (the most dominant) of the individuals involved in the repulsion-attraction process. Therefore, the social force, that is exerted between [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , is calculated by using the following new model: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the social force strength defined in (2) and [figure omitted; refer to PDF] is the unit vector pointing from [figure omitted; refer to PDF] to [figure omitted; refer to PDF] . Besides, [figure omitted; refer to PDF] is a randomly generated number between 1 and -1. Likewise, [figure omitted; refer to PDF] is the dominance function that calculates the dominance value of the most dominant individual from [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . In order to operate [figure omitted; refer to PDF] , all the individuals from [figure omitted; refer to PDF] are ranked according to their fitness values. The ranks are assigned so that the best individual receives the rank 0 (zero) whereas the worst individual obtains the rank [figure omitted; refer to PDF] . Therefore, the function [figure omitted; refer to PDF] is defined as follows: [figure omitted; refer to PDF] where the function rank( [figure omitted; refer to PDF] ) delivers the rank of the [figure omitted; refer to PDF] -individual. According to (8), [figure omitted; refer to PDF] yields a value within the interval [figure omitted; refer to PDF] . Its maximum value of one in [figure omitted; refer to PDF] is reached when either individual [figure omitted; refer to PDF] or [figure omitted; refer to PDF] is the best element of the population [figure omitted; refer to PDF] regarding their fitness values. On the other hand, a value close to zero is obtained when both individuals [figure omitted; refer to PDF] and [figure omitted; refer to PDF] possess quite bad fitness values. Figure 3 shows the behavior of [figure omitted; refer to PDF] considering 100 individuals. In the figure, it is assumed that [figure omitted; refer to PDF] represents one of the 99 individuals with ranks between 0 and 98 whereas [figure omitted; refer to PDF] is fixed to the element with the worst fitness value (rank 99).

Figure 3: Behavior of [figure omitted; refer to PDF] considering 100 individuals.

[figure omitted; refer to PDF]

Under the incorporation of [figure omitted; refer to PDF] in (7), social forces are magnified or weakened depending on the best fitness value (the most dominant) of the individuals involved in the repulsion-attraction process.

Finally, the total social force on each individual [figure omitted; refer to PDF] is modeled as the superposition of all of the pairwise interactions exerted over it: [figure omitted; refer to PDF] Therefore, the change of position [figure omitted; refer to PDF] is considered as the total social force experimented by [figure omitted; refer to PDF] as the superposition of all of the pairwise interactions. Thus, [figure omitted; refer to PDF] is defined as follows: [figure omitted; refer to PDF] After calculating the new positions [figure omitted; refer to PDF] of the population [figure omitted; refer to PDF] , the final positions [figure omitted; refer to PDF] must be calculated. The idea is to admit only the changes that guarantee an improvement in the search strategy. If the fitness value of [figure omitted; refer to PDF] is better than [figure omitted; refer to PDF] , then [figure omitted; refer to PDF] is accepted as the final solution. Otherwise, [figure omitted; refer to PDF] is retained. This procedure can be resumed by the following statement (considering a minimization problem): [figure omitted; refer to PDF] In order to illustrate the performance of the solitary operator, Figure 4 presents a simple example with the solitary operator being iteratively applied. A population of 50 different members ( [figure omitted; refer to PDF] ) is assumed which adopt a concentrated configuration as initial condition (Figure 4(a)). As a consequence of the social forces, the set of elements tends to distribute throughout the search space. Examples of different distributions are shown in Figures 4(b), 4(c), and 4(d) after applying 25, 50, and 100 different solitary operations, respectively.

Figure 4: Examples of different distributions. (a) Initial condition, (b) distribution after applying 25, (c) 50 and (d) 100 operations. The green asterisk represents the minimum value so far.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

3.2. Social Operation (B)

The social procedure represents the exploitation phase of the LS algorithm. Exploitation is the process of refining existent individuals within a small neighborhood in order to improve their solution quality.

The social procedure is a selective operation which is applied only to a subset [figure omitted; refer to PDF] of the final positions [figure omitted; refer to PDF] (where [figure omitted; refer to PDF] ). The operation starts by sorting [figure omitted; refer to PDF] with respect to fitness values, storing the elements in a temporary population [figure omitted; refer to PDF] . The elements in [figure omitted; refer to PDF] are sorted so that the best individual receives the position [figure omitted; refer to PDF] whereas the worst individual obtains the location [figure omitted; refer to PDF] . Therefore, the subset [figure omitted; refer to PDF] is integrated by only the first [figure omitted; refer to PDF] locations of [figure omitted; refer to PDF] (promising solutions). Under this operation, a subspace [figure omitted; refer to PDF] is created around each selected particle [figure omitted; refer to PDF] . The size of [figure omitted; refer to PDF] depends on the distance [figure omitted; refer to PDF] which is defined as follows: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the upper and lower bounds in the [figure omitted; refer to PDF] th dimension and [figure omitted; refer to PDF] is the number of dimensions of the optimization problem, whereas [figure omitted; refer to PDF] is a tuning factor. Therefore, the limits of [figure omitted; refer to PDF] can be modeled as follows: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the upper and lower bounds of the [figure omitted; refer to PDF] th dimension for the subspace [figure omitted; refer to PDF] , respectively.

Considering the subspace [figure omitted; refer to PDF] around each element [figure omitted; refer to PDF] , a set of [figure omitted; refer to PDF] new particles [figure omitted; refer to PDF] are randomly generated inside bounds fixed by (13). Once the [figure omitted; refer to PDF] samples are generated, the individual [figure omitted; refer to PDF] of the next population [figure omitted; refer to PDF] must be created. In order to calculate [figure omitted; refer to PDF] , the best particle [figure omitted; refer to PDF] , in terms of its fitness value from the [figure omitted; refer to PDF] samples (where [figure omitted; refer to PDF] ), is compared to [figure omitted; refer to PDF] . If [figure omitted; refer to PDF] is better than [figure omitted; refer to PDF] according to their fitness values, [figure omitted; refer to PDF] is updated with [figure omitted; refer to PDF] ; otherwise, [figure omitted; refer to PDF] is selected. The elements of [figure omitted; refer to PDF] that have not been processed by the procedure ( [figure omitted; refer to PDF] ) transfer their corresponding values to [figure omitted; refer to PDF] with no change.

The social operation is used to exploit only prominent solutions. According to the proposed method, inside each subspace [figure omitted; refer to PDF] , [figure omitted; refer to PDF] random samples are selected. Since the number of selected samples at each subspace is very small (typically [figure omitted; refer to PDF] ), the use of this operator substantially reduces the number of fitness function evaluations.

In order to demonstrate the social operation, a numerical example has been set by applying the proposed process to a simple function. Such function considers the interval of [figure omitted; refer to PDF] , [figure omitted; refer to PDF] whereas the function possesses one global maxima of value 8.1 at [figure omitted; refer to PDF] . Notice that [figure omitted; refer to PDF] and [figure omitted; refer to PDF] correspond to the axis coordinates (commonly [figure omitted; refer to PDF] and [figure omitted; refer to PDF] ). For this example, a final position population [figure omitted; refer to PDF] of six 2-dimensional members ( [figure omitted; refer to PDF] ) is assumed. Figure 5 shows the initial configuration of the proposed example, with the black points representing half of the particles with the best fitness values (the first three elements of [figure omitted; refer to PDF] , [figure omitted; refer to PDF] ) whereas the grey points [figure omitted; refer to PDF] correspond to the remaining individuals. From Figure 5, it can be seen that the social procedure is applied to all black particles ( [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] ) yielding two new random particles ( [figure omitted; refer to PDF] ), which are characterized by white points [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] , for each black point inside of their corresponding subspaces [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] . Considering the particle [figure omitted; refer to PDF] in Figure 7, the particle [figure omitted; refer to PDF] corresponds to the best particle ( [figure omitted; refer to PDF] ) from the two randomly generated particles (according to their fitness values) within [figure omitted; refer to PDF] . Therefore, the particle [figure omitted; refer to PDF] will substitute [figure omitted; refer to PDF] in the individual [figure omitted; refer to PDF] for the next generation, since it holds a better fitness value than [figure omitted; refer to PDF] .

Figure 5: Operation of the social procedure.

[figure omitted; refer to PDF]

The LS optimization procedure is defined over a bounded search space S . Search points that do not belong to such area are considered to be infeasible. However, during the evolution process, some candidate solutions could fall outside the search space. In the proposed approach, such infeasible solutions are arbitrarily placed with a random position inside the search space S .

3.3. Complete LS Algorithm

LS is a simple algorithm with only seven adjustable parameters: the strength of attraction [figure omitted; refer to PDF] , the attractive length [figure omitted; refer to PDF] , number of promising solutions [figure omitted; refer to PDF] , the population size [figure omitted; refer to PDF] , the tuning factor [figure omitted; refer to PDF] , the number of random samples [figure omitted; refer to PDF] , and the number of generations [figure omitted; refer to PDF] . The operation of LS is divided into three parts: initialization of the solitary and social operations. In the initialization ( [figure omitted; refer to PDF] ), the first population [figure omitted; refer to PDF] is produced. The values [figure omitted; refer to PDF] of each individual [figure omitted; refer to PDF] and each dimension [figure omitted; refer to PDF] are randomly and uniformly distributed between the prespecified lower initial parameter bound [figure omitted; refer to PDF] and the upper initial parameter bound [figure omitted; refer to PDF] : [figure omitted; refer to PDF] In the evolution process, the solitary (A) and social (B) operations are iteratively applied until the number of iterations [figure omitted; refer to PDF] has been reached. The complete LS procedure is illustrated in Algorithm 1.

Algorithm 1: Locust Search (LS) algorithm.

(1) Input : [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] and [figure omitted; refer to PDF] .

(2) Initialize [figure omitted; refer to PDF] ( [figure omitted; refer to PDF] = 0)

(3) until ( [figure omitted; refer to PDF] = [figure omitted; refer to PDF] )

(4) [figure omitted; refer to PDF] ( [figure omitted; refer to PDF] ) Solitary operator (Section 3.1)

(5) [figure omitted; refer to PDF] ( [figure omitted; refer to PDF] , [figure omitted; refer to PDF] ) Social operator (Section 3.2)

(6) [figure omitted; refer to PDF] = k + 1

(7) end until

3.4. Discussion about the LS Algorithm

Evolutionary algorithms (EA) have been widely employed for solving complex optimization problems. These methods are found to be more powerful than conventional methods that are based on formal logics or mathematical programming [46]. In the EA algorithm, search agents have to decide whether to explore unknown search positions or to exploit already tested positions in order to improve their solution quality. Pure exploration degrades the precision of the evolutionary process but increases its capacity to find new potentially solutions. On the other hand, pure exploitation allows refining existent solutions but adversely drives the process to local optimal solutions. Therefore, the ability of an EA to find a global optimal solution depends on its capacity to find a good balance between the exploitation of found-so-far elements and the exploration of the search space [47].

Most of swarm algorithms and other evolutionary algorithms tend to exclusively concentrate the individuals in the current best positions. Under such circumstances, such algorithms seriously limit their exploration-exploitation capacities.

Different to most of existent evolutionary algorithms, in the proposed approach, the modeled behavior explicitly avoids the concentration of individuals in the current best positions. Such fact allows not only to emulate the cooperative behavior of the locust colony in a good realistic way but also to incorporate a computational mechanism to avoid critical flaws that are commonly present in the popular evolutionary algorithms, such as the premature convergence and the incorrect exploration-exploitation balance.

It is important to emphasize that the proposed approach conducts two operators (solitary and social) within a single iteration. Such operators are similar to those that are used by other evolutionary methods such as ABC (employed bees, onlooker bees, and scout bees), AIS (clonal proliferation operator, affinity maturation operator, and clonal selection operator), and DE (mutation, crossover, and selection), which are all executed in a single iteration.

4. Numerical Experiments over Benchmark Functions

A comprehensive set of 13 functions, all collected from [48-50], has been used to test the performance of the LS approach as an optimization method. Tables 11 and 12 present the benchmark functions used in our experimental study. Such functions are classified into two different categories: unimodal test functions (Table 11) and multimodal test functions (Table 12). In these tables, [figure omitted; refer to PDF] is the function dimension; [figure omitted; refer to PDF] is the minimum value of the function, with [figure omitted; refer to PDF] being a subset of [figure omitted; refer to PDF] . The optimum locations ( [figure omitted; refer to PDF] ) for functions in Tables 11 and 12 are in [figure omitted; refer to PDF] , except for [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] with [figure omitted; refer to PDF] in [figure omitted; refer to PDF] and [figure omitted; refer to PDF] in [figure omitted; refer to PDF] . A detailed description of optimum locations is given in Tables 11 and 12.

We have applied the LS algorithm to 13 functions whose results have been compared to those produced by the Particle Swarm Optimization (PSO) method [27] and the Differential Evolution (DE) algorithm [28], both considered as the most popular algorithms for many optimization applications. In all comparisons, the population has been set to 40 ( [figure omitted; refer to PDF] ) individuals. The maximum iteration number for all functions has been set to 1000. Such stop criterion has been selected to maintain compatibility to similar works reported in the literature [48, 49].

The parameter setting for each of the algorithms in the comparison is described as follows:

(1) PSO: in the algorithm, [figure omitted; refer to PDF] while the inertia factor ( [figure omitted; refer to PDF] ) is decreased linearly from 0.9 to 0.2.

(2) DE: the DE/Rand/1 scheme is employed. The parameter settings follow the instructions in [28, 51]. The crossover probability is [figure omitted; refer to PDF] and the weighting factor is [figure omitted; refer to PDF] .

(3) In LS, [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are set to 0.6 and [figure omitted; refer to PDF] , respectively. Besides, [figure omitted; refer to PDF] is fixed to 20 ( [figure omitted; refer to PDF] /2), [figure omitted; refer to PDF] , [figure omitted; refer to PDF] whereas [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are configured to 1000 and 40, respectively. Once such parameters have been experimentally determined, they are kept for all experiments in this section.

In the comparison, three indexes are considered: the average best-so-far solution (ABS), the standard deviation (SD), and the number of function evaluations (NFE). The first two indexes assess the accuracy of the solution whereas the last one measures the computational cost. The average best-so-far solution (ABS) expresses the average value of the best function evaluations that have been obtained from 30 independent executions. The standard deviation (SD) indicates the dispersion of the ABS values. Evolutionary methods are, in general, complex pieces of software with several operators and stochastic branches. Therefore, it is difficult to conduct a complexity analysis from a deterministic perspective. Under such circumstances, it is more appropriate to use the number of function evaluations (NFE), just as it is used in the literature [52, 53], to evaluate and assess the computational effort (time) and the complexity among optimizers. It represents how many times an algorithm uses the objective function to evaluate the objective (fitness) function until the best solution of a determined execution has been found. Since the experiments require 30 different executions, the NFE index corresponds to the averaged value obtained from these executions. A small NFE value indicates that less time is needed to reach the global optimum.

4.1. Unimodal Test Functions

Functions [figure omitted; refer to PDF] to [figure omitted; refer to PDF] are unimodal functions. The results for unimodal functions over 30 runs are reported in Table 1 considering the average best-so-far solution (ABS), the standard deviation (SD), and the number of function evaluations (NFE). According to this table, LS provides better results than PSO and DE for all functions in terms of ABS and SD. In particular, the test yields the largest performance difference in functions [figure omitted; refer to PDF] . Such functions maintain a narrow curving valley that is hard to optimize, in case the search space cannot be explored properly and the direction changes cannot be kept up with [54]. For this reason, the performance differences are directly related to a better trade-off between exploration and exploitation that is produced by LS operators. In the practice, a main goal of an optimization algorithm is to find a solution as good as possible within a small time window. The computational cost for the optimizer is represented by its NFE values. According to Table 1, the NFE values that are obtained by the proposed method are smaller than its counterparts. Lower NFE values are more desirable since they correspond to less computational overload and, therefore, faster results. In the results perspective, it is clear that PSO and DE need more than 1000 generations in order to produce better solutions. However, this number of generations is considered in the experiments aiming for producing a visible contrast among the approaches. If the number of generations has been set to an exaggerated value, then all methods would converge to the best solution with no significant troubles.

Table 1: Minimization results from the benchmark functions test in Table 11 with [figure omitted; refer to PDF] . Maximum number of iterations = 1000.

		PSO	DE	LS
[figure omitted; refer to PDF]	ABS	1.66 × 10-1	6.27 × 10-3	4.55 × 10-4
SD	3.79 × 10-1	1.68 × 10-1	6.98 × 10-4
NFE	28,610	20,534	16,780
[figure omitted; refer to PDF]	ABS	4.83 × 10-1	2.02 × 10-1	5.41 × 10-3
SD	1.59 × 10-1	0.66	1.45 × 10-2
NFE	28,745	21,112	16,324
[figure omitted; refer to PDF]	ABS	2.75	5.72 × 10-1	1.61 × 10-3
SD	1.01	0.15	1.32 × 10-3
NFE	38,320	36,894	20,462
[figure omitted; refer to PDF]	ABS	1.84	0.11	1.05 × 10-2
SD	0.87	0.05	6.63 × 10-3
NFE	37,028	36,450	21,158
[figure omitted; refer to PDF]	ABS	3.07	2.39	4.11 × 10-2
SD	0.42	0.36	2.74 × 10-3
NFE	39,432	37,264	21,678
[figure omitted; refer to PDF]	ABS	6.36	6.51	5.88 × 10-2
SD	0.74	0.87	1.67 × 10-2
NFE	38,490	36,564	22,238
[figure omitted; refer to PDF]	ABS	6.14	0.12	2.71 × 10-2
SD	0.73	0.02	1.18 × 10-2
NFE	37,274	35,486	21,842

A nonparametric statistical significance proof known as Wilcoxon's rank sum test for independent samples [55, 56] has been conducted with an 5% significance level, over the "average best-so-far" data of Table 1. Table 2 reports the [figure omitted; refer to PDF] values produced by Wilcoxon's test for the pairwise comparison of the "average best-so-far" of two groups. Such groups are formed by LS versus PSO and LS versus DE. As a null hypothesis, it is assumed that there is no significant difference between mean values of the two algorithms. The alternative hypothesis considers a significant difference between the "average best-so-far" values of both approaches. All [figure omitted; refer to PDF] values reported in the table are less than 0.05 (5% significance level) which is a strong evidence against the null hypothesis, indicating that the LS results are statistically significant and that it has not occurred by coincidence (i.e., due to the normal noise contained in the process).

Table 2: [figure omitted; refer to PDF] values produced by Wilcoxon's test that compares LS versus PSO and DE over the "average best-so-far" values from Table 1.

LS versus	PSO	DE
[figure omitted; refer to PDF]	1.83 × 10-4	1.73 × 10-2
[figure omitted; refer to PDF]	3.85 × 10-3	1.83 × 10-4
[figure omitted; refer to PDF]	1.73 × 10-4	6.23 × 10-3
[figure omitted; refer to PDF]	2.57 × 10-4	5.21 × 10-3
[figure omitted; refer to PDF]	4.73 × 10-4	1.83 × 10-3
[figure omitted; refer to PDF]	6.39 × 10-5	2.15 × 10-3
[figure omitted; refer to PDF]	1.83 × 10-4	2.21 × 10-3

4.2. Multimodal Test Functions

Multimodal functions possess many local minima which make the optimization a difficult task to be accomplished. In multimodal functions, the results reflect the algorithm's ability to escape from local optima. We have applied the algorithms over functions [figure omitted; refer to PDF] to [figure omitted; refer to PDF] where the number of local minima increases exponentially as the dimension of the function increases. The dimension of such functions is set to 30. The results are averaged over 30 runs, with performance indexes being reported in Table 3 as follows: the average best-so-far solution (ABS), the standard deviation (SD), and the number of function evaluations (NFE). Likewise, [figure omitted; refer to PDF] values of the Wilcoxon signed-rank test of 30 independent runs are listed in Table 4. From the results, it is clear that LS yields better solutions than others algorithms for functions [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] , in terms of the indexes ABS and SD. However, for functions [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , LS produces similar results to DE. The Wilcoxon rank test results, that are presented in Table 6, confirm that LS performed better than PSO and DE considering four problems [figure omitted; refer to PDF] , whereas, from a statistical viewpoint, there is no difference between results from LS and DE for [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . According to Table 3, the NFE values obtained by the proposed method are smaller than those produced by other optimizers. The reason of this remarkable performance is associated with its two operators: (i) the solitary operator allows a better particle distribution in the search space, increasing the algorithm's ability to find the global optima and (ii) the use of the social operation provides a simple exploitation operator that intensifies the capacity of finding better solutions during the evolution process.

Table 3: Minimization results from the benchmark functions test in Table 12 with [figure omitted; refer to PDF] . Maximum number of iterations = 1000.

		PSO	DE	LS
[figure omitted; refer to PDF]	ABS	-6.7 × 103	-1.26 × 104	-1.26 × 104
SD	6.3 × 102	3.7 × 102	1.1 × 102
NFE	38,452	35,240	21,846
[figure omitted; refer to PDF]	ABS	14.8	4.01 × 10-1	2.49 × 10-3
SD	1.39	5.1 × 10-2	4.8 × 10-4
NFE	37,672	34,576	20,784
[figure omitted; refer to PDF]	ABS	14.7	4.66 × 10-2	2.15 × 10-3
SD	1.44	1.27 × 10-2	3.18 × 10-4
NFE	39,475	37,080	21,235
[figure omitted; refer to PDF]	ABS	12.01	1.15	1.47 × 10-4
SD	3.12	0.06	1.48 × 10-5
NFE	38,542	34,875	22,126
[figure omitted; refer to PDF]	ABS	6.87 × 10-1	3.74 × 10-1	5.58 × 10-3
SD	7.07 × 10-1	1.55 × 10-1	4.18 × 10-4
NFE	35,248	30,540	16,984
[figure omitted; refer to PDF]	ABS	1.87 × 10-1	1.81 × 10-2	1.78 × 10-2
SD	5.74 × 10-1	1.66 × 10-2	1.64 × 10-3
NFE	36,022	31,968	18,802

Table 4: [figure omitted; refer to PDF] values produced by Wilcoxon's test comparing LS versus PSO and DE over the "average best-so-far" values from Table 3.

LS versus	PSO	DE
[figure omitted; refer to PDF]	1.83 × 10-4	0.061
[figure omitted; refer to PDF]	1.17 × 10-4	2.41 × 10-4
[figure omitted; refer to PDF]	1.43 × 10-4	3.12 × 10-3
[figure omitted; refer to PDF]	6.25 × 10-4	1.14 × 10-3
[figure omitted; refer to PDF]	2.34 × 10-5	7.15 × 10-4
[figure omitted; refer to PDF]	4.73 × 10-4	0.071

5. Gaussian Mixture Modelling

In this section, the modeling of image histograms through Gaussian mixture models is presented. Let one consider an image holding [figure omitted; refer to PDF] gray levels [figure omitted; refer to PDF] whose distribution is defined by a histogram [figure omitted; refer to PDF] represented by the following formulation: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] denotes the number of pixels with gray level [figure omitted; refer to PDF] and Np the total number of pixels in the image. Under such circumstances, [figure omitted; refer to PDF] can be modeled by using a mixture [figure omitted; refer to PDF] of Gaussian functions of the form: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] symbolizes the number of Gaussian functions of the model whereas [figure omitted; refer to PDF] is the a priori probability of function [figure omitted; refer to PDF] . [figure omitted; refer to PDF] and [figure omitted; refer to PDF] represent the mean and standard deviation of the [figure omitted; refer to PDF] th Gaussian function, respectively. Furthermore, the constraint [figure omitted; refer to PDF] must be satisfied by the model. In order to evaluate the similarity between the image histogram and a candidate mixture model, the mean square error can be used as follows: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] represents the penalty associated with the constrain [figure omitted; refer to PDF] . Therefore, [figure omitted; refer to PDF] is considered as the objective function which must be minimized in the estimation problem. In order to illustrate the histogram modeling through a Gaussian mixture, Figure 6 presents an example, assuming three classes, that is, [figure omitted; refer to PDF] . Considering Figure 6(a) as the image histogram [figure omitted; refer to PDF] , the Gaussian mixture [figure omitted; refer to PDF] , that is shown in Figure 6(c), is produced by adding the Gaussian functions [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] in the configuration presented in Figure 6(b). Once the model parameters that better model the image histogram have been determined, the final step is to define the threshold values [figure omitted; refer to PDF] which can be calculated by simple standard methods, just as it is presented in [19-21].

Figure 6: Histogram approximation through a Gaussian mixture. (a) Original histogram, (b) configuration of the Gaussian components [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] , and (c) final Gaussian mixture [figure omitted; refer to PDF] .

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

Figure 7: Effect of the new objective function [figure omitted; refer to PDF] in the evaluation of Gaussian mixtures (candidate solutions). (a) Original histogram, (b) Gaussian mixture considering four classes, (c) penalization areas, and (d) Gaussian mixture of better quality solution.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

6. Segmentation Algorithm Based on LS

In the proposed method, the segmentation process is approached as an optimization problem. Computational optimization generally involves two distinct elements: (1) a search strategy to produce candidate solutions (individuals, particles, insects, locust, etc.) and (2) an objective function that evaluates the quality of each selected candidate solution. Several computational algorithms are available to perform the first element. The second element, where the objective function is designed, is unquestionably the most critical. In the objective function, it is expected to embody the full complexity of the performance, biases, and restrictions of the problem to be solved. In the segmentation problem, candidate solutions represent Gaussian mixtures. The objective function [figure omitted; refer to PDF] (17) is used as a fitness value to evaluate the similarity presented between the Gaussian mixture and the image histogram. Therefore, guided by the fitness values ( [figure omitted; refer to PDF] values), a set of encoded candidate solutions are evolved using the evolutionary operators until the best possible resemblance can be found.

Over the last decade, several algorithms based on evolutionary and swarm principles [19-22] have been proposed to solve the problem of segmentation through a Gaussian mixture model. Although these algorithms own certain good performance indexes, they present two important limitations. (1) They frequently obtain suboptimal approximations as a consequence of an inappropriate balance between exploration and exploitation in their search strategies. (2) They are based on the assumption that the number of Gaussians (classes) in the mixture is preknown and fixed; otherwise, they cannot work.

In order to eliminate such flaws, the proposed approach includes (A) a new search strategy and (B) the definition of a new objective function. For the search strategy, the LS method (Section 4) is adopted. Under LS, the concentration of individuals in the current best positions is explicitly avoided. Such fact allows reducing critical problems such as the premature convergence to suboptimal solutions and the incorrect exploration-exploitation balance.

6.1. New Objective Function [figure omitted; refer to PDF]

Previous segmentation algorithms based on evolutionary and swarm methods use (17) as objective function. Under these circumstances, each candidate solution (Gaussian mixture) is only evaluated in terms of its approximation with the image histogram.

Since the proposed approach aims to automatically select the number of Gaussian functions [figure omitted; refer to PDF] in the final mixture [figure omitted; refer to PDF] , the objective function must be modified. The new objective function [figure omitted; refer to PDF] is defined as follows: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is a scaling constant. The new objective function is divided into two parts. The first part [figure omitted; refer to PDF] evaluates the quality of each candidate solution in terms of its similarity with regard to the image histogram (17). The second part [figure omitted; refer to PDF] penalizes the overlapped area among Gaussian functions (classes), with [figure omitted; refer to PDF] being defined as follows: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] represent the number of classes and the gray levels, respectively. The parameters [figure omitted; refer to PDF] and [figure omitted; refer to PDF] symbolize the Gaussian functions [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , respectively, that are to be evaluated on the point [figure omitted; refer to PDF] (gray level) whereas the elements [figure omitted; refer to PDF] and [figure omitted; refer to PDF] represent the a priori probabilities of the Gaussian functions [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , respectively. Under such circumstances, mixtures with Gaussian functions that do not "positively" participate in the histogram approximation are severely penalized.

Figure 7 illustrates the effect of the new objective function [figure omitted; refer to PDF] in the evaluation of Gaussian mixtures (candidate solutions). From the image histogram (Figure 7(a)), it is evident that two Gaussian functions are enough to accurately approximate the original histogram. However, if the Gaussian mixture is modeled by using a greater number of functions (e.g., four as it is shown in Figure 7(b)), the original objective function [figure omitted; refer to PDF] is unable to obtain a reasonable approximation. Under the new objective function [figure omitted; refer to PDF] , the overlapped area among Gaussian functions (classes) is penalized. Such areas, in Figure 7(c), correspond to [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] , where [figure omitted; refer to PDF] represents the penalization value produced between the Gaussian function [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . Therefore, due to the penalization, the Gaussian mixture shown in Figures 7(b) and 7(c) provides a solution of low quality. On the other hand, the Gaussian mixture presented in Figure 7(d) maintains a low penalty; thus, it represents a solution of high quality. From Figure 7(d), it is easy to see that functions [figure omitted; refer to PDF] and [figure omitted; refer to PDF] can be removed from the final mixture. This elimination could be performed by a simple comparison with a threshold value [figure omitted; refer to PDF] , since [figure omitted; refer to PDF] and [figure omitted; refer to PDF] have a reduced amplitude ( [figure omitted; refer to PDF] ). Therefore, under [figure omitted; refer to PDF] , it is possible to find the reduced Gaussian mixture model starting from a considerable number of functions.

Since the proposed segmentation method is conceived as an optimization problem, the overall operation can be reduced to solve the formulation of (20) by using the LS algorithm: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] represent the probability, mean, and standard deviation of the class [figure omitted; refer to PDF] . It is important to remark that the new objective function [figure omitted; refer to PDF] allows the evaluation of a candidate solution in terms of the necessary number of Gaussian functions and its approximation quality. Under such circumstances, it can be used in combination with any other evolutionary method and not only with the proposed LS algorithm.

6.2. Complete Segmentation Algorithm

Once the new search strategy (LS) and objective function ( [figure omitted; refer to PDF] ) have been defined, the proposed segmentation algorithm can be summarized by Algorithm 2. The new algorithm combines operators defined by LS and operations for calculating the threshold values.

Algorithm 2: Segmentation LS algorithm.

(1) Input : [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] and [figure omitted; refer to PDF] .

(2) Initialize [figure omitted; refer to PDF] ( [figure omitted; refer to PDF] = 0)

(3) until ( [figure omitted; refer to PDF] = [figure omitted; refer to PDF] )

(4) [figure omitted; refer to PDF] ( [figure omitted; refer to PDF] ) Solitary operator (Section 3.1)

(5) [figure omitted; refer to PDF] ( [figure omitted; refer to PDF] , [figure omitted; refer to PDF] ) Social operator (Section 3.2)

(6) [figure omitted; refer to PDF] = k + 1

(7) end until

(8) Obtain [figure omitted; refer to PDF]

(9) Reduce [figure omitted; refer to PDF]

(10) Calculate the threshold values [figure omitted; refer to PDF] from the reduced model

(11) Use [figure omitted; refer to PDF] to segment the image

(Line 1) The algorithm sets the operative parameters [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] . They rule the behavior of the segmentation algorithm. (Line 2) Afterwards, the population [figure omitted; refer to PDF] is initialized considering [figure omitted; refer to PDF] different random Gaussian mixtures of [figure omitted; refer to PDF] functions. The idea is to generate an [figure omitted; refer to PDF] -random Gaussian mixture subjected to the constraints formulated in (20). The parameter [figure omitted; refer to PDF] must hold a high value in order to correctly segment all images (recall that the algorithm is able to reduce the Gaussian mixture to its minimum expression). (Line 3) Then, the Gaussian mixtures are evolved by using the LS operators and the new objective function [figure omitted; refer to PDF] . This process is repeated during [figure omitted; refer to PDF] cycles. (Line 8) After this procedure, the best Gaussian mixture [figure omitted; refer to PDF] is selected according to its objective function [figure omitted; refer to PDF] . (Line 9) Then, the Gaussian mixture [figure omitted; refer to PDF] is reduced by eliminating those functions whose amplitudes are lower than [figure omitted; refer to PDF] ( [figure omitted; refer to PDF] ). (Line 10) Then, the threshold values [figure omitted; refer to PDF] from the reduced model are calculated. (Line 11) Finally, the calculated [figure omitted; refer to PDF] values are employed to segment the image. Figure 8 shows a flowchart of the complete segmentation procedure.

Figure 8: Flowchart of the complete segmentation procedure.

[figure omitted; refer to PDF]

The proposed segmentation algorithm permits to automatically detect the number of segmentation partitions (classes). Furthermore, due to its remarkable search capacities, the LS method maintains a better accuracy than previous algorithms based on evolutionary principles. However, the proposed method presents two disadvantages: first, it is related to its implementation which in general is more complex than most of the other segmentators based on evolutionary basics. The second refers to the segmentation procedure of the proposed approach which does not consider any spatial pixel characteristics. As a consequence, pixels that may belong to a determined region due to its position are labeled as a part of another region due to its gray level intensity. Such a fact adversely affects the segmentation performance of the method.

7. Segmentation Results

This section analyses the performance of the proposed segmentation algorithm. The discussion is divided into three parts: the first one shows the performance of the proposed LS segmentation algorithm while the second presents a comparison between the proposed approach and others segmentation algorithms that are based on evolutionary and swam methods. The comparison mainly considers the capacities of each algorithm to accurately and robustly approximate the image histogram. Finally, the third part presents an objective evaluation of segmentation results produced by all algorithms that have been employed in the comparisons.

7.1. Performance of LS Algorithm in Image Segmentation

This section presents two experiments that analyze the performance of the proposed approach. Table 5 presents the algorithm's parameters that have been experimentally determined and kept for all the test images through all experiments.

Table 5: Parameter setup for the LS segmentation algorithm.

[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
0.6	0.2	20	1000	40	10	0.01	0.0001

Table 6: Results of the reduced Gaussian mixture for the first and the second image.

Parameters	First image	Second image
[figure omitted; refer to PDF]	0.004	0.032
[figure omitted; refer to PDF]	18.1	12.1
[figure omitted; refer to PDF]	8.2	2.9
[figure omitted; refer to PDF]	0.0035	0.0015
[figure omitted; refer to PDF]	69.9	45.1
[figure omitted; refer to PDF]	18.4	24.9
[figure omitted; refer to PDF]	0.01	0.006
[figure omitted; refer to PDF]	94.9	130.1
[figure omitted; refer to PDF]	8.9	17.9
[figure omitted; refer to PDF]	0.007	0.02
[figure omitted; refer to PDF]	163.1	167.2
[figure omitted; refer to PDF]	29.9	10.1

First Image. The first test considers the histogram shown by Figure 9(b) while Figure 9(a) presents the original image. After applying the proposed algorithm, just as it has been configured according to parameters in Table 5, a minimum model of four classes emerges after the start from Gaussian mixtures of 10 classes. Considering 30 independent executions, the averaged parameters of the resultant Gaussian mixture are presented in Table 6. One final Gaussian mixture (ten classes), which has been obtained by LS, is presented in Figure 10. Furthermore, the approximation of the reduced Gaussian mixture is also visually compared with the original histogram in Figure 10. On the other hand, Figure 11 presents the segmented image after calculating the threshold points.

Figure 9: (a) Original first image used on the first experiment and (b) its histogram.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

Figure 10: Gaussian mixture obtained by LS for the first image.

[figure omitted; refer to PDF]

Figure 11: Image segmented with the reduced Gaussian mixture.

[figure omitted; refer to PDF]

Second Image. For the second experiment, the image in Figure 12 is tested. The method aims to segment the image by using a reduced Gaussian mixture model obtained by the LS approach. After executing the algorithm according to parameters from Table 5, the resulting averaged parameters of the resultant Gaussian mixture are presented in Table 6. In order to assure consistency, the results are averaged considering 30 independent executions. Figure 13 shows the approximation quality that is obtained by the reduced Gaussian mixture model in (a) and the segmented image in (b).

Figure 12: (a) Original second image used on the first experiment and (b) its histogram.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

Figure 13: (a) Segmentation result obtained by LS for the first image and (b) the segmented image.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

7.2. Histogram Approximation Comparisons

This section discusses the comparison between the proposed segmentation algorithm and other evolutionary-segmentation methods that have been proposed in the literature. The set of methods used in the experiments includes the [figure omitted; refer to PDF] [19], [figure omitted; refer to PDF] [20], and [figure omitted; refer to PDF] [21]. These algorithms consider the combination between the original objective function [figure omitted; refer to PDF] (17) and an evolutionary technique such as Artificial Bee Colony (ABC), the Artificial Immune Systems (AIS), and the Differential Evolution (DE) [21], respectively. Since the proposed segmentation approach considers the combination of the new objective function [figure omitted; refer to PDF] (18) and the LS algorithm, it will be referred to as [figure omitted; refer to PDF] . The comparison focuses mainly on the capacities of each algorithm to accurately and robustly approximate the image histogram.

In the experiments, the populations have been set to 40 ( [figure omitted; refer to PDF] ) individuals. The maximum iteration number for all functions has been set to 1000. Such stop criterion has been considered to maintain compatibility to similar experiments that are reported in the literature [18]. The parameter setting for each of the segmentation algorithms in the comparison is described as follows:

(1) [figure omitted; refer to PDF] + ABC [19]: in the algorithm, its parameters are configured as follows: the abandonment [figure omitted; refer to PDF] , [figure omitted; refer to PDF] and [figure omitted; refer to PDF] .

(2) [figure omitted; refer to PDF] + AIS [20]: it presents the following parameters, [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] .

(3) [figure omitted; refer to PDF] + DE [21]: the DE/Rand/1 scheme is employed. The parameter settings follow the instructions in [21]. The crossover probability is [figure omitted; refer to PDF] and the weighting factor is [figure omitted; refer to PDF] .

(4) In [figure omitted; refer to PDF] , the method is set according to the values described in Table 5.

In order to conduct the experiments, a synthetic image is designed to be used as a reference in the comparisons. The main idea is to know in advance the exact number of classes (and their parameters) that are contained in the image so that the histogram can be considered as a ground truth. The synthetic image is divided into four sections. Each section corresponds to a different class which is produced by setting each gray level pixel [figure omitted; refer to PDF] to a value that is determined by the following model: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] represents the section, whereas [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the mean and the dispersion of the gray level pixels, respectively. The comparative study employs the image of [figure omitted; refer to PDF] that is shown in Figure 14(a) and the algorithm's parameters that have been presented in Table 7. Figure 14(b) illustrates the resultant histogram.

Table 7: Employed parameters for the design of the reference image.

	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
(1)	0.05	40	8
(2)	0.04	100	10
(3)	0.05	160	8
(4)	0.027	220	15

Figure 14: (a) Synthetic image used in the comparison and (b) its histogram.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

In the comparison, the discussion focuses on the following issues: first of all, accuracy; second, convergence; and third, computational cost.

Convergence. This section analyzes the approximation convergence when the number of classes that are used by the algorithm during the evolution process is different to the actual number of classes in the image. Recall that the proposed approach automatically finds the reduced Gaussian mixture which better adapts to the image histogram.

In the experiment, the methods, [figure omitted; refer to PDF] + ABC, [figure omitted; refer to PDF] + AIS, and [figure omitted; refer to PDF] + DE, are executed considering Gaussian mixtures composed of 8 functions. Under such circumstances, the number of classes to be detected is higher than the actual number of classes in the image. On the other hand, the proposed algorithm maintains the same configuration of Table 5. Therefore, it can detect and calculate up to ten classes ( [figure omitted; refer to PDF] ).

As a result, the techniques [figure omitted; refer to PDF] + ABC, [figure omitted; refer to PDF] + AIS, and [figure omitted; refer to PDF] + DE tend to overestimate the image histogram. This effect can be seen in Figure 15(a), where the resulting Gaussian functions are concentrated within actual classes. Such a behavior is a consequence of the evaluation that is considered by the original objective function [figure omitted; refer to PDF] , which privileges only the approximation between the Gaussian mixture and the image histogram. This effect can be graphically illustrated by Figure 16(a) that shows the pixel misclassification produced by the wrong segmentation of Figure 14(a). On the other hand, the proposed approach obtains a reduced Gaussian mixture model which allows the detection of each class from the actual histogram (see Figure 15(b)). As a consequence, the segmentation is significantly improved by eliminating the pixel misclassification, as it is shown by Figure 16(b).

Figure 15: Convergence results. (a) Convergence of the following methods: [figure omitted; refer to PDF] + ABC, [figure omitted; refer to PDF] + AIS, and [figure omitted; refer to PDF] + DE considering Gaussian mixtures of 8 classes. (b) Convergence of the proposed method (reduced Gaussian mixture).

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

Figure 16: Segmentation results obtained by (a) several methods including [figure omitted; refer to PDF] + ABC, [figure omitted; refer to PDF] + AIS, and [figure omitted; refer to PDF] + DE considering Gaussian mixtures of 8 classes and (b) the proposed method (reduced Gaussian mixture).

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

It is evident from Figure 15 that the techniques, [figure omitted; refer to PDF] + ABC, [figure omitted; refer to PDF] + AIS, and [figure omitted; refer to PDF] + DE, all need an a priori knowledge of the number of classes that are contained in the actual histogram in order to obtain a satisfactory result. On the other hand, the proposed algorithm is able to find a reduced Gaussian mixture whose classes coincide with the actual number of classes that are contained in the image histogram.

Accuracy. In this section, the comparison among the algorithms in terms of accuracy is reported. Most of the reported comparisons [19-26] are concerned about comparing the parameters of the resultant Gaussian mixtures by using real images. Under such circumstances, it is difficult to consider a clear reference in order to define a meaningful error. Therefore, the image defined in Figure 14 has been used in the experiments because its construction parameters are clearly defined in Table 7.

Since the parameter values from Table 7 act as ground truth, a simple averaged difference between them and the values that are computed by each algorithm could be used as comparison error. However, as each parameter maintains different active intervals, it is necessary to express the differences in terms of percentage. Therefore, if [figure omitted; refer to PDF] is the parametric difference and [figure omitted; refer to PDF] the ground truth parameter, the percentage error [figure omitted; refer to PDF] can be defined as follows: [figure omitted; refer to PDF] In the segmentation problem, each Gaussian mixture represents a [figure omitted; refer to PDF] -dimensional model where each dimension corresponds to a Gaussian function of the optimization problem to be solved. Since each Gaussian function possesses three parameters [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] , the complete number of parameters is [figure omitted; refer to PDF] dimensions. Therefore, the final error [figure omitted; refer to PDF] produced by the final Gaussian mixture is [figure omitted; refer to PDF] where [figure omitted; refer to PDF] .

In order to compare accuracy, the algorithms, [figure omitted; refer to PDF] + ABC, [figure omitted; refer to PDF] + AIS, [figure omitted; refer to PDF] + DE, and the proposed approach are all executed over the image shown by Figure 14(a). The experiment aims to estimate the Gaussian mixture that better approximates the actual image histogram. Methods [figure omitted; refer to PDF] + ABC, [figure omitted; refer to PDF] + AIS, and [figure omitted; refer to PDF] + DE consider Gaussian mixtures composed of 4 functions ( [figure omitted; refer to PDF] ). In case of the [figure omitted; refer to PDF] method, although the algorithm finds a reduced Gaussian mixture of four functions, it is initially set with ten functions ( [figure omitted; refer to PDF] ). Table 8 presents the final Gaussian mixture parameters and the final error [figure omitted; refer to PDF] . The final Gaussian mixture parameters have been averaged over 30 independent executions in order to assure consistency. A close inspection of Table 8 reveals that the proposed method is able to achieve the smallest error ( [figure omitted; refer to PDF] ) in comparison to the other algorithms. Figure 16 presents the histogram approximations that are produced by each algorithm whereas Figure 17 shows their correspondent segmented images. Both illustrations present the median case obtained throughout 30 runs. Figure 18 exhibits that [figure omitted; refer to PDF] + ABC, [figure omitted; refer to PDF] + AIS, and [figure omitted; refer to PDF] + DE present different levels of misclassifications which are nearly absent in the proposed approach case.

Table 8: Results of the reduced Gaussian mixture in terms of accuracy.

Algorithm	Gaussian function	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
J + ABC	(1)	0.052	44.5	6.4	11.79%
(2)	0.084	98.12	12.87
(3)	0.058	163.50	8.94
(4)	0.025	218.84	17.5
J + AIS	(1)	0.072	31.01	6.14	22.01%
(2)	0.054	88.52	12.21
(3)	0.039	149.21	9.14
(4)	0.034	248.41	13.84
J + DE	(1)	0.041	35.74	7.86	13.57%
(2)	0.036	90.57	11.97
(3)	0.059	148.47	9.01
(4)	0.020	201.34	13.02
[figure omitted; refer to PDF] + LS	(1)	0.049	40.12	7.5	3.98%
(2)	0.041	102.04	10.4
(3)	0.052	168.66	8.3
(4)	0.025	110.92	15.2

Figure 17: Approximation results in terms of accuracy. (a) [figure omitted; refer to PDF] + ABC, (b) [figure omitted; refer to PDF] + AIS, (c) [figure omitted; refer to PDF] + DE, and (d) the proposed [figure omitted; refer to PDF] approach.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

Figure 18: Segmentation results in terms of accuracy. (a) [figure omitted; refer to PDF] + ABC, (b) [figure omitted; refer to PDF] + AIS, (c) [figure omitted; refer to PDF] + DE, and (d) the proposed [figure omitted; refer to PDF] approach.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

Computational Cost. The experiment aims to measure the complexity and the computing time spent by the [figure omitted; refer to PDF] + ABC, the [figure omitted; refer to PDF] + AIS, the [figure omitted; refer to PDF] + DE, and the [figure omitted; refer to PDF] algorithm while calculating the parameters of the Gaussian mixture in benchmark images (see Figures 19(a)-19(d)). [figure omitted; refer to PDF] + ABC, [figure omitted; refer to PDF] + AIS, and [figure omitted; refer to PDF] + DE consider Gaussian mixtures that are composed of 4 functions ( [figure omitted; refer to PDF] ). In case of the [figure omitted; refer to PDF] method, although the algorithm finds a reduced Gaussian mixture of four functions despite being initialized with ten functions ( [figure omitted; refer to PDF] ), Table 9 shows the averaged measurements after 30 experiments. It is evident that the [figure omitted; refer to PDF] + ABC and [figure omitted; refer to PDF] + DE are the slowest to converge (iterations) and the [figure omitted; refer to PDF] + AIS shows the highest computational cost (time elapsed) because it requires operators which demand long times. On the other hand, the [figure omitted; refer to PDF] shows an acceptable trade-off between its convergence time and its computational cost. Therefore, although the implementation of [figure omitted; refer to PDF] in general requires more code than most of other evolution-based segmentators, such a fact is not reflected in the execution time. Finally, Figure 19 below shows the segmented images as they are generated by each algorithm. It can be seen that the proposed approach generate more homogeneous regions whereas [figure omitted; refer to PDF] + ABC, [figure omitted; refer to PDF] + AIS, and [figure omitted; refer to PDF] + DE present several artifacts that are produced by an incorrect pixel classification.

Table 9: Iterations and time requirements of the [figure omitted; refer to PDF] + ABC, the [figure omitted; refer to PDF] + AIS, the [figure omitted; refer to PDF] + DE. and the [figure omitted; refer to PDF] + LS algorithm as they are applied to segment benchmark images (see Figure 17).

Iterations	(a)	(b)	(c)	(d)
Time elapsed
[figure omitted; refer to PDF] + ABC	855	833	870	997
2.72 s	2.70 s	2.73 s	3.1 s
[figure omitted; refer to PDF] + AIS	725	704	754	812
1.78 s	1.61 s	1.41 s	2.01 s
[figure omitted; refer to PDF] + DE	657	627	694	742
1.25 s	1.12 s	1.45 s	1.88 s
[figure omitted; refer to PDF] + LS	314	298	307	402
0.98 s	0.84 s	0.72 s	1.02 s

Figure 19: Images employed in the computational cost analysis.

[figure omitted; refer to PDF]

7.3. Performance Evaluation of the Segmentation Results

This section presents an objective evaluation of segmentation results that are produced by all algorithms in the comparisons. The ill-defined nature of the segmentation problem makes the evaluation of a candidate algorithm difficult [57]. Traditionally, the evaluation has been conducted by using some supervised criteria [58] which are based on the computation of a dissimilarity measure between a segmentation result and a ground truth image. Recently, the use of unsupervised measures has substituted supervised indexes for the objective evaluation of segmentation results [59]. They enable the quantification of the quality of a segmentation result without a priori knowledge (ground truth image).

Evaluation Criteria. In this paper, the unsupervised index ROS proposed by Chabrier et al. [60] has been used to objectively evaluate the performance of each candidate algorithm. This index evaluates the segmentation quality in terms of the homogeneity within segmented regions and the heterogeneity among the different regions. ROS can be computed as follows: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] quantifies the homogeneity within segmented regions. Similarly, [figure omitted; refer to PDF] measures the disparity among the regions. A segmentation result [figure omitted; refer to PDF] is considered better than [figure omitted; refer to PDF] , if [figure omitted; refer to PDF] . The interregion homogeneity characterized by [figure omitted; refer to PDF] is calculated considering the following formulation: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] represents the number of partitions in which the image has been segmented. [figure omitted; refer to PDF] symbolizes the number of pixels contained in the partition [figure omitted; refer to PDF] whereas [figure omitted; refer to PDF] considers the number of pixels that integrate the complete image. Similarly, [figure omitted; refer to PDF] represents the standard deviation from the partition [figure omitted; refer to PDF] . On the other hand, disparity among the regions [figure omitted; refer to PDF] is computed as follows: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the average gray level in the partition [figure omitted; refer to PDF] .

Experimental Protocol. In the comparison of segmentation results, a set of four classical images has been chosen to integrate the experimental set (Figure 20). The segmentation methods used in the comparison are [figure omitted; refer to PDF] + ABC [19], [figure omitted; refer to PDF] + AIS [20], and [figure omitted; refer to PDF] + DE [21].

Figure 20: Experimental set used in the evaluation of the segmentation results.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

From all segmentation methods used in the comparison, the proposed [figure omitted; refer to PDF] algorithm is the only one that has the capacity to automatically detect the number of segmentation partitions (classes). In order to conduct a fair comparison, all algorithms have been proved by using the same number of partitions. Therefore, in the experiments, the [figure omitted; refer to PDF] segmentation algorithm is firstly applied to detect the best possible number of partitions [figure omitted; refer to PDF] . Once we obtained the number of partitions [figure omitted; refer to PDF] , the rest of the algorithms were configured to approximate the image histogram with this number of classes.

Figure 21 presents the segmentation results obtained by each algorithm considering the experimental set from Figure 20. On the other hand, Table 10 shows the evaluation of the segmentation results in terms of the ROS index. Such values represent the averaged measurements after 30 executions. From them, it can be seen that the proposed [figure omitted; refer to PDF] method obtains the best ROS indexes. Such values indicate that the proposed algorithm maintains the best balance between the homogeneity within segmented regions and the heterogeneity among the different regions. From Figure 21, it can be seen that the proposed approach generates more homogeneous regions whereas [figure omitted; refer to PDF] + ABC, [figure omitted; refer to PDF] + AIS, and [figure omitted; refer to PDF] + DE present several artifacts that are produced by an incorrect pixel classification.

Table 10: Evaluation of the segmentation results in terms of the ROS index.

Number of classesImage	[figure omitted; refer to PDF] (a)	[figure omitted; refer to PDF] (b)	[figure omitted; refer to PDF] (c)	[figure omitted; refer to PDF] (d)
[figure omitted; refer to PDF] + ABC	0.534	0.328	0.411	0.457
[figure omitted; refer to PDF] + AIS	0.522	0.321	0.427	0.437
[figure omitted; refer to PDF] + DE	0.512	0.312	0.408	0.418
[figure omitted; refer to PDF] + LS	0.674	0.401	0.514	0.527

Table 11: Unimodal test functions.

Test function	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	0
[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	0
[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	0
[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	0
[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	0
[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	0
[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	0

Table 12: Multimodal test functions.

Test function	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	0
[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	0
[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	0
[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	0
[figure omitted; refer to PDF] [figure omitted; refer to PDF]
[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	0

Figure 21: Segmentation results using in the evaluation.

[figure omitted; refer to PDF]

8. Conclusions

Despite the fact that several evolutionary methods have been successfully applied to image segmentation with interesting results, most of them have exhibited two important limitations: (1) they frequently obtain suboptimal results (misclassifications) as a consequence of an inappropriate balance between exploration and exploitation in their search strategies; (2) the number of classes is fixed and known in advance.

In this paper, a new swarm algorithm for the automatic image segmentation, called the Locust Search (LS), has been presented. The proposed method eliminates the typical flaws presented by previous evolutionary approaches by combining a novel evolutionary method with the definition of a new objective function that appropriately evaluates the segmentation quality with respect to the number of classes. In order to illustrate the proficiency and robustness of the proposed approach, several numerical experiments have been conducted. Such experiments have been divided into two parts. First, the proposed LS method has been compared to other well-known evolutionary techniques on a set of benchmark functions. In the second part, the performance of the proposed segmentation algorithm has been compared to other segmentation methods based on evolutionary principles. The results in both cases validate the efficiency of the proposed technique with regard to accuracy and robustness.

Several research directions will be considered for future work such as the inclusion of other indexes to evaluate similarity between a candidate solution and the image histogram, the consideration of spatial pixel characteristics in the objective function, the modification of the evolutionary LS operators to control the exploration-exploitation balance, and the conversion of the segmentation procedure into a multiobjective problem.

Acknowledgment

The present research has been supported under the Grant CONACYT CB 181053.

Conflict of Interests

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

References

[1] H. Zhang, J. E. Fritts, S. A. Goldman, "Image segmentation evaluation: a survey of unsupervised methods," Computer Vision and Image Understanding , vol. 110, no. 2, pp. 260-280, 2008.

[2] T. Uemura, G. Koutaki, K. Uchimura, "Image segmentation based on edge detection using boundary code," International Journal of Innovative Computing , vol. 7, pp. 6073-6083, 2011.

[3] L. Wang, H. Wu, C. Pan, "Region-based image segmentation with local signed difference energy," Pattern Recognition Letters , vol. 34, no. 6, pp. 637-645, 2013.

[4] N. Otsu, "A thresholding selection method from gray-level histogram," IEEE Transactions on Systems, Man and Cybernetics , vol. 9, no. 1, pp. 62-66, 1979.

[5] R. Peng, P. K. Varshney, "On performance limits of image segmentation algorithms," Computer Vision and Image Understanding , vol. 132, pp. 24-38, 2015.

[6] M. A. Balafar, "Gaussian mixture model based segmentation methods for brain MRI images," Artificial Intelligence Review , vol. 41, no. 3, pp. 429-439, 2014.

[7] G. J. McLachlan, S. Rathnayake, "On the number of components in a Gaussian mixture model," Data Mining and Knowledge Discovery , vol. 4, no. 5, pp. 341-355, 2014.

[8] D. Oliva, V. Osuna-Enciso, E. Cuevas, G. Pajares, M. Perez-Cisneros, D. Zaldivar, "Improving segmentation velocity using an evolutionary method," Expert Systems with Applications , vol. 42, no. 14, pp. 5874-5886, 2015.

[9] Z.-W. Ye, M.-W. Wang, W. Liu, S.-B. Chen, "Fuzzy entropy based optimal thresholding using bat algorithm," Applied Soft Computing , vol. 31, pp. 381-395, 2015.

[10] S. Sarkar, S. Das, S. Sinha Chaudhuri, "A multilevel color image thresholding scheme based on minimum cross entropy and differential evolution," Pattern Recognition Letters , vol. 54, pp. 27-35, 2015.

[11] A. K. Bhandari, A. Kumar, G. K. Singh, "Modified artificial bee colony based computationally efficient multilevel thresholding for satellite image segmentation using Kapur's, Otsu and Tsallis functions," Expert Systems with Applications , vol. 42, no. 3, pp. 1573-1601, 2015.

[12] H. Permuter, J. Francos, I. Jermyn, "A study of Gaussian mixture models of color and texture features for image classification and segmentation," Pattern Recognition , vol. 39, no. 4, pp. 695-706, 2006.

[13] A. P. Dempster, A. P. Laird, D. B. Rubin, "Maximum likelihood from incomplete data via the EM algorithm," Journal of the Royal Statistical Society Series B , vol. 39, no. 1, pp. 1-38, 1977.

[14] Z. Zhang, C. Chen, J. Sun, K. L. Chan, "EM algorithms for Gaussian mixtures with split-and-merge operation," Pattern Recognition , vol. 36, no. 9, pp. 1973-1983, 2003.

[15] H. Park, S.-I. Amari, K. Fukumizu, "Adaptive natural gradient learning algorithms for various stochastic models," Neural Networks , vol. 13, no. 7, pp. 755-764, 2000.

[16] H. Park, T. Ozeki, "Singularity and slow convergence of the em algorithm for gaussian mixtures," Neural Processing Letters , vol. 29, no. 1, pp. 45-59, 2009.

[17] L. Gupta, T. Sortrakul, "A Gaussian-mixture-based image segmentation algorithm," Pattern Recognition , vol. 31, no. 3, pp. 315-325, 1998.

[18] V. Osuna-Enciso, E. Cuevas, H. Sossa, "A comparison of nature inspired algorithms for multi-threshold image segmentation," Expert Systems with Applications , vol. 40, no. 4, pp. 1213-1219, 2013.

[19] E. Cuevas, F. Sencion, D. Zaldivar, M. Perez-Cisneros, H. Sossa, "A multi-threshold segmentation approach based on artificial bee colony optimization," Applied Intelligence , vol. 37, no. 3, pp. 321-336, 2012.

[20] E. Cuevas, V. Osuna-Enciso, D. Zaldivar, M. Perez-Cisneros, H. Sossa, "Multithreshold segmentation based on artificial immune systems," Mathematical Problems in Engineering , vol. 2012, 2012.

[21] E. Cuevas, D. Zaldivar, M. Perez-Cisneros, "A novel multi-threshold segmentation approach based on differential evolution optimization," Expert Systems with Applications , vol. 37, no. 7, pp. 5265-5271, 2010.

[22] D. Oliva, E. Cuevas, G. Pajares, D. Zaldivar, V. Osuna, "A multilevel thresholding algorithm using electromagnetism optimization," Neurocomputing , vol. 139, pp. 357-381, 2014.

[23] D. Oliva, E. Cuevas, G. Pajares, D. Zaldivar, M. Perez-Cisneros, "Multilevel thresholding segmentation based on harmony search optimization," Journal of Applied Mathematics , vol. 2013, 2013.

[24] E. Cuevas, D. Zaldivar, M. Perez-Cisneros, "Seeking multi-thresholds for image segmentation with Learning Automata," Machine Vision and Applications , vol. 22, no. 5, pp. 805-818, 2011.

[25] K. C. Tan, S. C. Chiam, A. A. Mamun, C. K. Goh, "Balancing exploration and exploitation with adaptive variation for evolutionary multi-objective optimization," European Journal of Operational Research , vol. 197, no. 2, pp. 701-713, 2009.

[26] G. Chen, C. P. Low, Z. Yang, "Preserving and exploiting genetic diversity in evolutionary programming algorithms," IEEE Transactions on Evolutionary Computation , vol. 13, no. 3, pp. 661-673, 2009.

[27] J. Kennedy, R. Eberhart, "Particle swarm optimization," in Proceedings of the IEEE International Conference on Neural Networks, pp. 1942-1948, December 1995.

[28] R. Storn, K. Price, "Differential Evolution-a simple and efficient adaptive scheme for global optimisation over continuous spaces,", no. TR-95-012, ICSI, Berkeley, Calif, USA, 1995.

[29] Y. Wang, B. Li, T. Weise, J. Wang, B. Yuan, Q. Tian, "Self-adaptive learning based particle swarm optimization," Information Sciences , vol. 181, no. 20, pp. 4515-4538, 2011.

[30] J. Tvrdik, "Adaptation in differential evolution: a numerical comparison," Applied Soft Computing Journal , vol. 9, no. 3, pp. 1149-1155, 2009.

[31] H. Wang, H. Sun, C. Li, S. Rahnamayan, J.-s. Pan, "Diversity enhanced particle swarm optimization with neighborhood search," Information Sciences , vol. 223, pp. 119-135, 2013.

[32] W. Gong, Á. Fialho, Z. Cai, H. Li, "Adaptive strategy selection in differential evolution for numerical optimization: an empirical study," Information Sciences , vol. 181, no. 24, pp. 5364-5386, 2011.

[33] D. M. Gordon, "The organization of work in social insect colonies," Complexity , vol. 8, no. 1, pp. 43-46, 2002.

[34] S. Kizaki, M. Katori, "Stochastic lattice model for locust outbreak," Physica A: Statistical Mechanics and its Applications , vol. 266, no. 1-4, pp. 339-342, 1999.

[35] S. M. Rogers, D. A. Cullen, M. L. Anstey, M. Burrows, E. Despland, T. Dodgson, T. Matheson, S. R. Ott, K. Stettin, G. A. Sword, S. J. Simpson, "Rapid behavioural gregarization in the desert locust, Schistocerca gregaria entails synchronous changes in both activity and attraction to conspecifics," Journal of Insect Physiology , vol. 65, pp. 9-26, 2014.

[36] C. M. Topaz, A. J. Bernoff, S. Logan, W. Toolson, "A model for rolling swarms of locusts," European Physical Journal: Special Topics , vol. 157, no. 1, pp. 93-109, 2008.

[37] C. M. Topaz, M. R. D'Orsogna, L. Edelstein-Keshet, A. J. Bernoff, "Locust dynamics: behavioral phase change and swarming," PLoS Computational Biology , vol. 8, no. 8, 2012.

[38] G. Oster, E. Wilson Caste and Ecology in the Social Insects , Princeton University Press, Princeton, NJ, USA, 1978.

[39] B. Hölldobler, E. O. Wilson Journey to the Ants: A Story of Scientific Exploration , 1994.

[40] B. Hölldobler, E. O. Wilson The Ants , Harvard University Press, Cambridge, Mass, USA, 1990.

[41] S. Tanaka, Y. Nishide, "Behavioral phase shift in nymphs of the desert locust, Schistocerca gregaria : special attention to attraction/avoidance behaviors and the role of serotonin," Journal of Insect Physiology , vol. 59, no. 1, pp. 101-112, 2013.

[42] E. Gaten, S. J. Huston, H. B. Dowse, T. Matheson, "Solitary and gregarious locusts differ in circadian rhythmicity of a visual output neuron," Journal of Biological Rhythms , vol. 27, no. 3, pp. 196-205, 2012.

[43] I. Benaragama, J. R. Gray, "Responses of a pair of flying locusts to lateral looming visual stimuli," Journal of Comparative Physiology A , vol. 200, no. 8, pp. 723-738, 2014.

[44] M. G. Sergeev, "Distribution patterns of grasshoppers and their Kin in the boreal zone," Psyche , vol. 2011, 2011.

[45] S. O. Ely, P. G. N. Njagi, M. O. Bashir, S. El-Tom El-Amin, A. Hassanali, "Diel behavioral activity patterns in adult solitarious desert locust, Schistocerca gregaria (Forskål)," Psyche , vol. 2011, 2011.

[46] X.-S. Yang Nature-Inspired Metaheuristic Algorithms , Luniver Press, Beckington, UK, 2008.

[47] E. Cuevas, A. Echavarria, M. A. Ramirez-Ortegon, "An optimization algorithm inspired by the States of Matter that improves the balance between exploration and exploitation," Applied Intelligence , vol. 40, no. 2, pp. 256-272, 2014.

[48] M. M. Ali, C. Khompatraporn, Z. B. Zabinsky, "A numerical evaluation of several stochastic algorithms on selected continuous global optimization test problems," Journal of Global Optimization , vol. 31, no. 4, pp. 635-672, 2005.

[49] R. Chelouah, P. Siarry, "Continuous genetic algorithm designed for the global optimization of multimodal functions," Journal of Heuristics , vol. 6, no. 2, pp. 191-213, 2000.

[50] F. Herrera, M. Lozano, A. M. Sanchez, "A taxonomy for the crossover operator for real-coded genetic algorithms: an experimental study," International Journal of Intelligent Systems , vol. 18, no. 3, pp. 309-338, 2003.

[51] E. Cuevas, D. Zaldivar, M. Perez-Cisneros, M. Ramirez-Ortegon, "Circle detection using discrete differential evolution optimization," Pattern Analysis & Applications , vol. 14, no. 1, pp. 93-107, 2011.

[52] A. Sadollah, H. Eskandar, A. Bahreininejad, J. H. Kim, "Water cycle algorithm with evaporation rate for solving constrained and unconstrained optimization problems," Applied Soft Computing , vol. 30, pp. 58-71, 2015.

[53] M. S. Kiran, H. Hakli, M. Gunduz, H. Uguz, "Artificial bee colony algorithm with variable search strategy for continuous optimization," Information Sciences , vol. 300, pp. 140-157, 2015.

[54] F. Kang, J. Li, Z. Ma, "Rosenbrock artificial bee colony algorithm for accurate global optimization of numerical functions," Information Sciences , vol. 181, no. 16, pp. 3508-3531, 2011.

[55] F. Wilcoxon, "Individual comparisons by ranking methods," Biometrics Bulletin , vol. 1, no. 6, pp. 80-83, 1945.

[56] S. Garcia, D. Molina, M. Lozano, F. Herrera, "A study on the use of non-parametric tests for analyzing the evolutionary algorithms' behaviour: a case study on the CEC'2005 Special Session on Real Parameter Optimization," Journal of Heuristics , vol. 15, no. 6, pp. 617-644, 2009.

[57] R. Unnikrishnan, C. Pantofaru, M. Hebert, "Toward objective evaluation of image segmentation algorithms," IEEE Transactions on Pattern Analysis and Machine Intelligence , vol. 29, no. 6, pp. 929-944, 2007.

[58] R. Unnikrishnan, C. Pantofaru, M. Hebert, "A measure for objective evaluation of image segmentation algorithms," in Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR '05), pp. 34, San Diego, Calif, USA, June 2005.

[59] Y. J. Zhang, "A survey on evaluation methods for image segmentation," Pattern Recognition , vol. 29, no. 8, pp. 1335-1346, 1996.

[60] S. Chabrier, B. Emile, C. Rosenberger, H. Laurent, "Unsupervised performance evaluation of image segmentation," EURASIP Journal on Advances in Signal Processing , vol. 2006, 2006.

Appendix

Appendix

List of Benchmark Functions

In Table 11, [figure omitted; refer to PDF] is the dimension of function, [figure omitted; refer to PDF] is the minimum value of the function, and [figure omitted; refer to PDF] is a subset of [figure omitted; refer to PDF] . The optimum location ( [figure omitted; refer to PDF] ) for functions in Table 11 is in [figure omitted; refer to PDF] , except for [figure omitted; refer to PDF] with [figure omitted; refer to PDF] in [figure omitted; refer to PDF] .

The optimum locations ( [figure omitted; refer to PDF] ) for functions in Table 12 are in [figure omitted; refer to PDF] , except for [figure omitted; refer to PDF] in [figure omitted; refer to PDF] and [figure omitted; refer to PDF] - [figure omitted; refer to PDF] in [figure omitted; refer to PDF] .