# A Multilevel Image Thresholding Method Based on Subspace Elimination Optimization

Authors: Jun Qin; ChuTing Wang; GuiHe Qin
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

Multilevel thresholding is to find the thresholds to segment the image with grey levels. Usually, the thresholds are so determined that some indicator functions of the segmented image are optimized. To improve the computational efficiency, we presented an optimization method for multilevel thresholding. First, the solution space is divided into subspaces. Second, the subspaces are searched to obtain their current local optimal value. Third, the subspaces that are of worse current optimal value are eliminated. Then, the next round of elimination is exerted in the remainder. The elimination is repeated until only one subspace is left and its optimal value is taken as the global optimum. In principle, any random search algorithm can be used to find the local optimum in a subspace block because it is a strategy to enhance the searching efficiency through eliminating hopeless regions as early as possible, rather than to improve the searching algorithm itself. To verify its performance, taking PSO (Particle swarm optimization) as the basic searching algorithm of subspaces, the presented method is applied to Otsu’s and Kapur’s multilevel thresholding of four different kinds of digital images. The presented method is compared with PSO, and it behaves better in efficiency.


---

## Body

## 1. Introduction

Image segmentation is one of the most important techniques in image processing and is an essential procedure in almost all image and video applications. Image segmentation is to classify the pixels of an image into different sets by such attribute as grey level. In an optimum segmentation, pixels should be so divided that some objective function value is optimized which are defined according to significance of the segmentation or the application area of the image. One kind of the most used segmentation techniques is image thresholding based on the image grey-level histogram [1–4]. Otsu [5] and Kapur [6] optimal thresholding segmentation are typical image grey-level histogram methods which take the variance and entropy between the classes as objective function, respectively.The k thresholds optimum segmentation is to find a group of grey level (t1,t2,⋯tk) to make some index optimized. It can be considered as an optimization problem to find the best objective function (the index) point in a k-dimension space in which all dimensions are the grey levels. So, when the grey-level histogram methods as Otsu’s and Kapur’s methods are extended to multilevel thresholding problems, its efficiency becomes very low because the computational complexity of this kind method increases exponentially due to its exhaustive search [7, 8]. The computational complex is one of the major reasons that so much different heuristic optimal techniques have been used in multilevel thresholding [9–11]. Generally, heuristic algorithms are inspired by such things as natural phenomena, physical laws, biological social activities, and evolutionary processes [10–12]. Many multilevel thresholding methods have been developed from these kinds of algorithms, and almost every new heuristic is used in multilevel thresholding not long after it is proposed as Human Mental Search (HMS) [13–15]. These optimization methods combined with Otsu’s and Kapur’s segmentation index are widely used in multithresholding segmentation [16–21].Among the heuristic methods of multilevel thresholding, many are based on particle swam optimization (PSO) algorithms for its simple and practical features as well as computational efficiency [11, 17, 21, 22]. PSO algorithm has been proven to be a powerful competitor to other heuristic algorithms in such optimization application as multilevel thresholding [23] and many enhanced PSO-based multilevel thresholding methods have been proposed.Most of the PSO-based multilevel thresholding methods directly use basic or improved PSO as optimization algorithm to search thresholds. Huang presented a method for multilevel thresholding segmentation based on QPSO (quantum particle swarms algorithm) algorithm in which the thresholds of Otsu’s method are founded with the QPSO algorithm [24]. The speed and accuracy are improved compared to traditional methods. Chen proposed a hybrid algorithm based on a self-adaptive PSO to optimize the thresholds of Ostu’s method [9]. Tang improved the parameters and evolutional process of basic PSO and used it in multilevel image thresholding based on maximum entropy [25]. The proposed method could get better optimal thresholds and more ideal segmentation result with less computation cost and had good search stability compared with the thresholding based on basic PSO in the experiments. In [7], Zhang improved particle swarm optimization with the improvements of particle’s best fitness value as the inertia weight of basic PSO and used the PSO to improve the selecting speed of the thresholds of Otsu’s method. Frank Jiang developed a particle swarm optimization with wavelet mutation named (PSOWM, particle swarm optimization with wavelet mutation) and used it in the multilevel thresholding image segmentation [26]. Their basic idea was to optimize the thresholds for the image segmentation to maximize the total entropy of the segmented image. According to their experiment, the new algorithm led to better thresholds, produced more accurate segmentation results for images with multiple attributes, and well balanced search stability and computational cost. In [8], Liu presented PSO improvement mechanisms which were named adaptive inertia (AI) and adaptive population (AP) to modify the basic PSO and used the modified PSO in multilevel thresholding to increase search efficiency and convergence speed. Adithya Alva proposed a Tsallis entropy-based multilevel thresholding technique using an image enhancing technique and an improvised algorithm called Half-Life Constant Particle Swarm Optimization (HCPSO) [27]. They tested the algorithm on different standard test images and professed that objective values and PSNR values are improved compared with the basic PSO and other optimization algorithms such as Genetic Algorithm (GA) and Bacterial Foraging (BF) algorithm. In [28], Habba Maryam proposed a multilevel thresholding method with multiobjective particle swarm optimization (MOPSO) in which the best thresholds determination process was seen as a multiobjective optimization problem. The objective functions of the MOPSO are based on the entropy calculation of the segmented regions. The proposed method was tested on different images (filtered by a low pass filter) and the experimental results showed its efficiency.In field of multilevel thresholding, there are also methods based on local particle swarm optimization (LPSO). Lin proposed a PSO-based multilevel thresholding method which embedded a scheme based on maximum entropy and uniformity into the particle swarm optimization with subswarm structure to find the optimal thresholds [29]. The swarm was divided into several subswarms to get local optimal results with an objective function based on maximal entropy and uniformity. In the method, just one-swarm particle was used to replace k-swarm (k is the number of thresholds) particles in order to upgrade the computation performance. The local optimal solutions were used to update global parameters in the global swarm. Through iterations updating, the velocities and locations of particles, the near optimal thresholds could be calculated out based on the objective functions of maximum entropy and uniformity. The hybrid algorithms which combine several heuristics or combine heuristics with mathematical programming are also widely used in multilevel thresholding of image segmentation [30–32]. These hybrid algorithms usually could improve the performance in solving large complex optimization problems [33].According to “no free lunch (NFL) theorems,” no optimization algorithm is the best for all problems. The method of optimization should take full consideration of the characteristics of the problem to solve [34–36]. Therefore, when designing an optimization method for multilevel thresholding, we should consider its characteristics. Among the various improved PSO algorithms for threshold segmentation, the improvement is aimed at the algorithm itself with little consideration of the specificity of multithresholding. Threshold segmentation is an exhaustive optimization problem in which the dimension of decision variables space, usually equal to the number of the thresholds, is not very high. And, for the actual images, especially that are composed of foreground background, some parts of the decision variable space are not likely to have the optimal thresholds.To improve efficiency of the exploring work, when searching for a target, such as prospecting, people usually partition the whole area to subregions, eliminate the subregions that the target exists in them with small possibility after a preliminary surveying, and refine the exploration of the remainder regions. Inspired by the idea as well as the characteristics of multilevel thresholding, we presented a subspace elimination optimization (SSEO) method and applied it in multilevel thresholding. First of all, the decision variables space is divided into subspace blocks. Second, the subspace blocks are searched to obtain their local optimal value estimation, respectively. Third, the subspace blocks with poor estimated optimal value are eliminated. Then, the next round of elimination is exerted in the winners of last round. The elimination is repeated until only one block is remained. The optimal value of it is taken as the global one. Taking PSO (Particle swarm optimization) as the basic searching method of subspace blocks, the presented method is applied in Otsu’s and Kapur’s multilevel thresholding of four different kinds of digital images. The experiment results show the effectiveness in enhancement of efficiency of heuristic optimal methods such as PSO. This method is to optimize the search space rather than to enhance the search method itself. So, any heuristic algorithms can be used to search subspace blocks. In optimal value exploration, eliminating hopeless parts of the solution space as early as possible and paying exploitation in the hopeful parts are beneficial to improve the efficiency of optimization process.


## 2. Subspace Elimination Optimization Method

### 2.1. Idea and Analysis

The canonical PSO algorithm searches the whole decision space with a swarm to find the optimal objective function value [37]. Extensive studies reveal that PSO’s performance mainly depends upon its two characteristics: exploration and exploitation [38, 39]. The process could be construed as that every particle tries to find the further best position according to its own historical optimal value and the one having achieved by the whole swarm. It is based on the conviction that the real optimal point could be found in the positions around the current optimal points with a large probability and should be explored with more dense exploitation. The basic steps are as follows.(a) The particle (searcher) constantly explores the objective function value of the neighbourhood points in the solution space and remembers the best location searched so far.(b) The particles communicate with others to learn the global best position so far.(c) All the particles move toward the global best position from their own local best position.Step (a) implies that the direct basic obtainment method of the knowledge on the function is for every particle to explore their local area. Step (b) gives the mechanism for the particles to know global knowledge with sharing local knowledge through communications. Step (c) is based on “the conviction that the global optimal point could be found in the positions around the current optimal points with a large probability.” For a specific optimization problem, the characteristics of the objective function determine this probability. (c) is also a major factor for PSO to fall into the local extremum [40]. Many improved PSO methods to avoid falling into the local optimum are based on adjustment of the fly speed calculating strategy with the expense of increasing computational complexity or decreasing accuracy.Local version PSO (LPSO) is a kind of these enhanced algorithms. Some researches manifest that LPSO is adopted in exploration, while GPSO is beneficial for exploitation [41]. Some studies have shown that LPSO is beneficial to avoid falling into local optimum, but the computational complexity is relatively higher [42, 43]. But, after analysing a large number of optimization problem samples, AP Engelbrecht deems that LPSO and GPSO (global version PSO) perform very similarly with respect to solution accuracy and LPSO performs slightly better than GPSO with respect to consistency [44]. Their main point is that “if the objective is to find the most optimal algorithm for a specific optimization problem, then the neighbourhood topology used has to be included as one of the parameters tuned. However, when the objective is to determine if a specific change to PSO provides improved performance, then it is arbitrary whether GPSO or LPSO is used” [44]. This is consistent with “no free lunch” (NFL) theorems mentioned above.In resource exploration or competition, the elimination method is often used to find the best objective. The area of exploration is partitioned into some parts or groups first. Then, all parts are explored and the poor parts are eliminated. And the potential parts are to be explored and eliminated further. The course is constantly executed until the objective emerges. Based on the observation, when searching the optimum, we can divide the solution space into some subspaces, and constantly eliminating the poor parts through exploring their performance with a heuristic algorithm as PSO. The remainders enter next-round elimination competition. Round by round, when only one subspace is left, the best position in the subspace is the optimum. The exploiting intensity for each subspace increases with the rounds. Many studies indicate that PSO with a sparse neighbour topology yields favourable results on complicated multimodal problems, while PSO with a dense neighbour topology offers promising solutions on simple unimodal problems [45]. The presented idea means that the space is explored with sparse swarm at the beginning, and denser swarm is used when less subspaces are left where the objective function behave simpler in modal for the space is much smaller.In principle, any heuristic algorithm can be used in searching the optimum in a subspace. With no loss of generality, we search the subspace with PSO and compare the presented method with the canonical PSO. This method is essentially not an improved PSO method. PSO is used independently in every subspace. There is no interconnection between the PSO algorithms in subspace. The essential strategy is exploration and elimination of subspaces. PSO is employed in the paper based on the reason that it is suitable for the method because it is simple and effective as well as its performance could be adjusted expediently by the swarm size and search times. The method is used in multilevel thresholding with objective function of Otus’s maximum between class variance and Kapur’s maximum entropy being objective functions, respectively, to illustrate its performance. The experiment results expounded that the method could improve the search efficiency compared with the method used in subspace.


### 2.2. Subspace Division and Identification

For a function defined in n-dimensional spacefx-, how to find the optimum of fx- is its optimization problem. In a digital quantity (such as grey level of digital image) space, decision variable x- can be defined as:(1)x-∈x0,x1,…xn-1(2)xj∈0,1,2,…Mi-1,j=0,1,2,…n-1Mj is a positive integer. For example, when xj express the grey level of a pixel of digital image with 8-bit binary number, Mj=255.First of all, solution space is divided into some subspaces which have the uniform size and same dimension as that of the solution space. To divide the solution space, each dimension is evenly divided into several segments, and the segment length injth dimension is marked as Sj. Each subspace is identified as shown in Table 1 with vector whose components are the segment order number starting from the origin point. A subspace block is labelled as B(k0,k1,⋯kn-1), in which kj is the segment order number in jth dimension.Table 1
Subspace division and identification.
Subspace label Regionin dimension 0 Regionin dimension 1 ⋯ Regionin dimension n-1 B ( 0 , 0 , ⋯ 0 ) [0,S0) [0,S1) ⋯ [0,Sn-1)  B ( 1 , 0 , ⋯ 0 ) [1∙S0, 2∙S0) [0,S1) ⋯ [0,Sn-1)  B ( 2 , 0 , ⋯ 0 ) [2∙S0, 3∙S0) [0,S1) ⋯ [0,Sn-1)  ⋯ ⋯ ⋯ ⋯ ⋯  B ( N 0 - 1 , 0 , ⋯ 0 ) [(N0-1)∙S0, 3∙S0) [0,S1) ⋯ [0,Sn-)  ⋯ ⋯ ⋯ ⋯ ⋯  B ( k 0 , k 1 , ⋯ k n - 1 ) [k0∙S0,(k0+1)∙S0) [k1∙S1,(k1+1)∙S1) ⋯ [kn-1∙Sn-1,(kn-1+1)∙Sn-1)  ⋯ ⋯ ⋯ ⋯ ⋯  B ( N 0 - 1 , N 1 - 1 , ⋯ N n - 1 - 1 ) [(N0-1)∙S0, N0∙S0) [(N1-1)∙S1, N1∙S1) ⋯ [(Nn-1-1)∙Sn-1, Nn-1∙Sn-1)The region ofBk0,k1,⋯kn-1 in jth dimension is(3)kj∙Sj,kj+1∙Sj,j=0,1,2,…n-1Figure 1 shows that the examples of dimension n are, respectively, 1 and 2, and every dimension is partitioned into 4 segments.Figure 1
Examples of subspace division and identification.
(a)
One dimension with 4 segments (b)
Two dimensions with 4 segments


### 2.3. Exploration of a Subspace

In every subspace, any optimization method can be adopted in principle. Here, PSO is used. Its concrete velocity calculation formula is defined as in equation (5):(4)Xijt+1=kj∙Sj,Xijt+Vijt+1<kj∙SjXijt+Vijt+1,kj∙Sj<Xijt+Vijt+1<kj+1∙Sjkj+1∙Sj,Xijt+Vijt+1>kj+1∙Sj(5)Vijt+1=ω∙Vijt+c1∙r1t∙Pijt-Xijt+c2∙r2t∙Gjt-Xijtwhere xijt+1 is component of the position of particle i in dimension j at time t+1, Pijt, Gjt are, respectively, the components of the optimal position of the particle i and the global optimal position of the subspace block in dimension j at time t, c1 and c2 are coefficients, r1t and r2t are random number in [0,1], and i=1,2,⋯NS, j=0,2,⋯n-1.


### 2.4. The Algorithm

According to the above analysis and the subspace division and identification, the subspace elimination optimization (SSEO) algorithm can be deduced as follows.Step 1.
Set initial subspace search particle numberNS, search times CS, search round number t=0, and the state of all subspace blocks to “active.”Step 2.
Evenly distribute particles in every subspaceB(k0,k1,⋯kn-1) as(6)Pi,jBk0,k1,⋯kn-1t←kj∙Sj+1NS+1∙Sjwhere Pi,jB(k0,k1,⋯kn-1)(t) is jth dimension component of the position of the particle i att.Step 3.
Calculate the initial optimum value ofB(k0,k1,⋯kn-1):(7)fmaxBk0,k1,⋯kn-1=maxi=1,2,⋯,NS⁡fPiBk0,k1,⋯kn-1tStep 4.
Determine the initial best position inB(k0,k1,⋯kn-1):(8)GBk0,k1,⋯kn-1t=ARCfmaxBk0,k1,⋯kn-1Step 5.
Search all the “active” subspace blocksCS times with PSO by NS particles.Step 6.
Assign the new optimum value ofB(k0,k1,⋯kn-1) to fmaxB(k0,k1,⋯kn-1).Step 7.
Record the coordinate of the new best position inB(k0,k1,⋯kn-1):(9)GBk0,k1,⋯kn-1t=ARCfmaxBk0,k1,⋯kn-1Step 8.
Sort all the “active” blocks according to their optimum.Step 9.
Eliminate the inferior blocks (marked as “passive”), and if only one block remains at the state “active”, then go to Step11.Step 10.
Adjust the subspace block search parametersNS and CS; go to Step 5 to execute next-round elimination.Step 11.
Take the optimum value of the last “active” subspace block as the global one and its coordinate as optimal point.


## 2.1. Idea and Analysis

The canonical PSO algorithm searches the whole decision space with a swarm to find the optimal objective function value [37]. Extensive studies reveal that PSO’s performance mainly depends upon its two characteristics: exploration and exploitation [38, 39]. The process could be construed as that every particle tries to find the further best position according to its own historical optimal value and the one having achieved by the whole swarm. It is based on the conviction that the real optimal point could be found in the positions around the current optimal points with a large probability and should be explored with more dense exploitation. The basic steps are as follows.(a) The particle (searcher) constantly explores the objective function value of the neighbourhood points in the solution space and remembers the best location searched so far.(b) The particles communicate with others to learn the global best position so far.(c) All the particles move toward the global best position from their own local best position.Step (a) implies that the direct basic obtainment method of the knowledge on the function is for every particle to explore their local area. Step (b) gives the mechanism for the particles to know global knowledge with sharing local knowledge through communications. Step (c) is based on “the conviction that the global optimal point could be found in the positions around the current optimal points with a large probability.” For a specific optimization problem, the characteristics of the objective function determine this probability. (c) is also a major factor for PSO to fall into the local extremum [40]. Many improved PSO methods to avoid falling into the local optimum are based on adjustment of the fly speed calculating strategy with the expense of increasing computational complexity or decreasing accuracy.Local version PSO (LPSO) is a kind of these enhanced algorithms. Some researches manifest that LPSO is adopted in exploration, while GPSO is beneficial for exploitation [41]. Some studies have shown that LPSO is beneficial to avoid falling into local optimum, but the computational complexity is relatively higher [42, 43]. But, after analysing a large number of optimization problem samples, AP Engelbrecht deems that LPSO and GPSO (global version PSO) perform very similarly with respect to solution accuracy and LPSO performs slightly better than GPSO with respect to consistency [44]. Their main point is that “if the objective is to find the most optimal algorithm for a specific optimization problem, then the neighbourhood topology used has to be included as one of the parameters tuned. However, when the objective is to determine if a specific change to PSO provides improved performance, then it is arbitrary whether GPSO or LPSO is used” [44]. This is consistent with “no free lunch” (NFL) theorems mentioned above.In resource exploration or competition, the elimination method is often used to find the best objective. The area of exploration is partitioned into some parts or groups first. Then, all parts are explored and the poor parts are eliminated. And the potential parts are to be explored and eliminated further. The course is constantly executed until the objective emerges. Based on the observation, when searching the optimum, we can divide the solution space into some subspaces, and constantly eliminating the poor parts through exploring their performance with a heuristic algorithm as PSO. The remainders enter next-round elimination competition. Round by round, when only one subspace is left, the best position in the subspace is the optimum. The exploiting intensity for each subspace increases with the rounds. Many studies indicate that PSO with a sparse neighbour topology yields favourable results on complicated multimodal problems, while PSO with a dense neighbour topology offers promising solutions on simple unimodal problems [45]. The presented idea means that the space is explored with sparse swarm at the beginning, and denser swarm is used when less subspaces are left where the objective function behave simpler in modal for the space is much smaller.In principle, any heuristic algorithm can be used in searching the optimum in a subspace. With no loss of generality, we search the subspace with PSO and compare the presented method with the canonical PSO. This method is essentially not an improved PSO method. PSO is used independently in every subspace. There is no interconnection between the PSO algorithms in subspace. The essential strategy is exploration and elimination of subspaces. PSO is employed in the paper based on the reason that it is suitable for the method because it is simple and effective as well as its performance could be adjusted expediently by the swarm size and search times. The method is used in multilevel thresholding with objective function of Otus’s maximum between class variance and Kapur’s maximum entropy being objective functions, respectively, to illustrate its performance. The experiment results expounded that the method could improve the search efficiency compared with the method used in subspace.


## 2.2. Subspace Division and Identification

For a function defined in n-dimensional spacefx-, how to find the optimum of fx- is its optimization problem. In a digital quantity (such as grey level of digital image) space, decision variable x- can be defined as:(1)x-∈x0,x1,…xn-1(2)xj∈0,1,2,…Mi-1,j=0,1,2,…n-1Mj is a positive integer. For example, when xj express the grey level of a pixel of digital image with 8-bit binary number, Mj=255.First of all, solution space is divided into some subspaces which have the uniform size and same dimension as that of the solution space. To divide the solution space, each dimension is evenly divided into several segments, and the segment length injth dimension is marked as Sj. Each subspace is identified as shown in Table 1 with vector whose components are the segment order number starting from the origin point. A subspace block is labelled as B(k0,k1,⋯kn-1), in which kj is the segment order number in jth dimension.Table 1
Subspace division and identification.
Subspace label Regionin dimension 0 Regionin dimension 1 ⋯ Regionin dimension n-1 B ( 0 , 0 , ⋯ 0 ) [0,S0) [0,S1) ⋯ [0,Sn-1)  B ( 1 , 0 , ⋯ 0 ) [1∙S0, 2∙S0) [0,S1) ⋯ [0,Sn-1)  B ( 2 , 0 , ⋯ 0 ) [2∙S0, 3∙S0) [0,S1) ⋯ [0,Sn-1)  ⋯ ⋯ ⋯ ⋯ ⋯  B ( N 0 - 1 , 0 , ⋯ 0 ) [(N0-1)∙S0, 3∙S0) [0,S1) ⋯ [0,Sn-)  ⋯ ⋯ ⋯ ⋯ ⋯  B ( k 0 , k 1 , ⋯ k n - 1 ) [k0∙S0,(k0+1)∙S0) [k1∙S1,(k1+1)∙S1) ⋯ [kn-1∙Sn-1,(kn-1+1)∙Sn-1)  ⋯ ⋯ ⋯ ⋯ ⋯  B ( N 0 - 1 , N 1 - 1 , ⋯ N n - 1 - 1 ) [(N0-1)∙S0, N0∙S0) [(N1-1)∙S1, N1∙S1) ⋯ [(Nn-1-1)∙Sn-1, Nn-1∙Sn-1)The region ofBk0,k1,⋯kn-1 in jth dimension is(3)kj∙Sj,kj+1∙Sj,j=0,1,2,…n-1Figure 1 shows that the examples of dimension n are, respectively, 1 and 2, and every dimension is partitioned into 4 segments.Figure 1
Examples of subspace division and identification.
(a)
One dimension with 4 segments (b)
Two dimensions with 4 segments


## 2.3. Exploration of a Subspace

In every subspace, any optimization method can be adopted in principle. Here, PSO is used. Its concrete velocity calculation formula is defined as in equation (5):(4)Xijt+1=kj∙Sj,Xijt+Vijt+1<kj∙SjXijt+Vijt+1,kj∙Sj<Xijt+Vijt+1<kj+1∙Sjkj+1∙Sj,Xijt+Vijt+1>kj+1∙Sj(5)Vijt+1=ω∙Vijt+c1∙r1t∙Pijt-Xijt+c2∙r2t∙Gjt-Xijtwhere xijt+1 is component of the position of particle i in dimension j at time t+1, Pijt, Gjt are, respectively, the components of the optimal position of the particle i and the global optimal position of the subspace block in dimension j at time t, c1 and c2 are coefficients, r1t and r2t are random number in [0,1], and i=1,2,⋯NS, j=0,2,⋯n-1.


## 2.4. The Algorithm

According to the above analysis and the subspace division and identification, the subspace elimination optimization (SSEO) algorithm can be deduced as follows.Step 1.
Set initial subspace search particle numberNS, search times CS, search round number t=0, and the state of all subspace blocks to “active.”Step 2.
Evenly distribute particles in every subspaceB(k0,k1,⋯kn-1) as(6)Pi,jBk0,k1,⋯kn-1t←kj∙Sj+1NS+1∙Sjwhere Pi,jB(k0,k1,⋯kn-1)(t) is jth dimension component of the position of the particle i att.Step 3.
Calculate the initial optimum value ofB(k0,k1,⋯kn-1):(7)fmaxBk0,k1,⋯kn-1=maxi=1,2,⋯,NS⁡fPiBk0,k1,⋯kn-1tStep 4.
Determine the initial best position inB(k0,k1,⋯kn-1):(8)GBk0,k1,⋯kn-1t=ARCfmaxBk0,k1,⋯kn-1Step 5.
Search all the “active” subspace blocksCS times with PSO by NS particles.Step 6.
Assign the new optimum value ofB(k0,k1,⋯kn-1) to fmaxB(k0,k1,⋯kn-1).Step 7.
Record the coordinate of the new best position inB(k0,k1,⋯kn-1):(9)GBk0,k1,⋯kn-1t=ARCfmaxBk0,k1,⋯kn-1Step 8.
Sort all the “active” blocks according to their optimum.Step 9.
Eliminate the inferior blocks (marked as “passive”), and if only one block remains at the state “active”, then go to Step11.Step 10.
Adjust the subspace block search parametersNS and CS; go to Step 5 to execute next-round elimination.Step 11.
Take the optimum value of the last “active” subspace block as the global one and its coordinate as optimal point.


## 3. Multilevel Thresholding with the Algorithm

Threshold segmentation can be understood as an optimization problem in the space made up of grey level. That is, to find the levels to class the pixels by which some objective functions achieve the optimum. A typical classification of threshold segmentation is based on objective function, such as minimum error, maximum interclass variance, maximum entropy, etc. Generally, Otsu’s method, Kapur’s method, Tsallis entropy, and minimum cross entropy are the best methods in thresholding based on optimizing the objective function [3]. Otsu’s method and Kapur’s method are, respectively, the most typical examples of the methods based on variance and entropy. Based on the objective functions of the two methods, we achieve multilevel thresholding with the proposed algorithm.K thresholds digital image segmentation could correspond to an optimal problem in k dimension space:(10)0,L-1x0,L-1⋯x0,L-1where L is the grey level.The coordinate components of the optimum point are the thresholds.For example, 3 thresholds in Otsu maximum variance (σ2) segmentation for an 8-bit digital image correspond to the optimization to find the maximum σ2 in the space 0,255x0,255x0,255. If the coordinate of the optimum point is (t1,t2,t3), the three thresholds are t1,t2,t3, respectively.For an L grey-level digital image G which is ofNC column and NL row, its pixel number is NLxNC, and grey-level set is {0,1,2,L-1}. The probability with a pixel having a grey level gi is(11)Pi=giNLxNC,i∈0,1,2,⋯L-1where pixels in image G are divided into k+1 classes C0, C1, …, and Ck according to their grey level with k thresholds t1,t2⋯tk. The grey level range of class C0 is 0,t1-1, C1 is t1,t2-1,…, and Ck is tk,L-1.The probabilityωi of a pixel belonging to the class Ci and the grey level probability mean μi of class Ci can be calculated as follows:(12)ωi=∑j=titi+1-1Pj,μi=∑j=titi+1-1j∙Pjωiwhere i=0,1,2,…,k, t0=0, and tk+1=L.The grey level mean value of the image is(13)μ=∑i=0L-1i∙PiThe Otsu interclass variance [5] is(14)σ2=∑i=0kωi∙μi-μ2The Kapur entropy value [6] is(15)Ψ=-∑i=0t1-1piω0ln⁡piω0-∑i=t1t2-1piω1ln⁡piω1-⋯∑i=tkL-1piωkln⁡piωkThe Otsu maximum between class variance multilevel thresholding is to find the thresholds t1,t2⋯tk with which σ2 is maximized, and the Kapur maximum entropy multilevel thresholding is to find the thresholds t1,t2⋯tk with which Ψ is maximized.


## 4. Results and Discussion

To test the effectiveness of the proposed method, Otus’s and Kapur’s multilevel thresholding of four images are conducted with the SSEO method. The basic optimal algorithm used in subspaces is PSO, and the results of the method are also compared with that when PSO is directly applied to calculate the thresholds.The four original images are from the USC-SIPI image library as illustrated in Figure2. The experiment is carried out on a PC with 64-bit 2.70GHz CPU (Intel(R) Core(TM)i5-6400 CPU@2.70GHz), 8.00GB main memory, and Windows10. Software platform through which the algorithm is programed is MATLAB R2016a.Figure 2
Original images.
(a)
Female (256∗256) (b)
Jelly beans (256∗256) (c)
San Diego (512∗512) (d)
Sailboat (512∗512)Taking the Otsu maximum interclass variance and Kapur maximum entropy as the objective functions, 1 to 4 thresholds of the four mages are calculated with exhaustion method, PSO, and SSEO method, respectively. The parametersc1=2, c2=2 in the PSO algorithm. The results of exhaustion method are taken as the reference standard. The end criterion of PSO is that the searched results are close enough to them. In 1 threshold situation, the difference is not bigger than 1, 2 threshold situation is 2, and so on. Therefore, this is the fastest situation that to obtain the result for PSO algorithm. The algorithm of searching the subspace is the same as the global PSO with a little revise as (4) and much less particles. End criterion of SSEO is search times. Table 2 shows the partition of solution space. Tables 3 and 4 show the threshold computation results with the three methods for Otsu maximum interclass variance and Kapur maximum entropy segmentation. Accounting for the stochastic characteristics, the results of PSO and SSEO are average values of 10 times calculation results.Table 2
Partition of solution space.
Number of thresholds (dimension) 1 2 3 4 Partition of solution space 5 3 ∗ 3 = 9 4 ∗ 4 ∗ 4 = 64 4 ∗ 4 ∗ 4 ∗ 4 = 256Table 3
Otsu’s maximum between class variance multilevel thresholding.
Image Number of thresholds Index Exhaustion method PSO SSEO Female 1 Variance 1145.60 1145.25 1145.60 Thresholds 75 75 74 Time(s) 0.085460 0.108066 0.122840 Variance of the Time  7.58789E-06 8.52760E-06 2 Variance 1347.10 1346.51 1346.83 Thresholds 58 110 59 110 58 110 Time(s) 0.280200 0.138602 0.155361 Variance of the Time  1.11086E-04 2.80965E-06 3 Variance 1444.70 1443.30 1444.04 Thresholds 35 69 116 35 69 115 34 69 115 Time(s) 18.351619 0.243077 0.227970 Variance of the Time  2.55672E-03 2.30173E-05 4 Variance 1487.40 1486.08 1486.85 Thresholds 33 63 99 138 33 64 100 140 33 63 99 139 Time(s) 1132.261655 0.766454 0.623988 Variance of the Time  1.37653E-01 1.66319E-05  Jelly beans 1 Variance 1236.50 1236.26 1236.50 Thresholds 139 139 138 Time(s) 0.091220 0.109814 0.121777 Variance of the Time  9.40414E-06 5.22426E-06 2 Variance 1349.60 1349.42 1349.37 Thresholds 89 155 89 155 103 166 Time(s) 0.291574 0.142287 0.155299 Variance of the Time  2.38194E-04 2.01919E-06 3 Variance 1399.80 1399.14 1399.21 Thresholds 85 145 186 84 145 185 85 145 185 Time(s) 18.154217 7.573480 0.228355 Variance of the Time  2.40161E+02 7.18563E-06 4 Variance 1428.90 1428.18 1428.25 Thresholds 76 119 160 188 74 119 161 187 76 119 159 188 Time(s) 1244.274622 0.902986 0.620868 Variance of the Time  2.61322E-01 4.99843E-05  San Diego 1 Variance 1712.50 1712.22 1712.44 Thresholds 141 141 140 Time(s) 0.096458 0.118310 0.130571 Variance of the Time  7.54120E-06 1.14140E-05 2 Variance 2074.00 2073.39 2073.53 Thresholds 106 165 105 165 105 164 Time(s) 0.300369 0.161998 0.168375 Variance of the Time  4.54634E-04 1.30487E-06 3 Variance 2188.80 2187.52 2187.48 Thresholds 90 128 176 91 127 175 90 128 175 Time(s) 18.114168 0.270145 0.245001 Variance of the Time  1.74912E-02 1.38766E-05 4 Variance 2249.20 2247.51 2248.10 Thresholds 76 108 143 184 76 108 143 183 75 108 144 184 Time(s) 1109.358181 0.899595 0.646769 Variance of the Time  4.54533E-01 1.28512E-05  Sailboat 1 Variance 3689.50 3689.42 3689.48 Thresholds 126 125 125 Time(s) 0.097266 0.114995 0.132179 Variance of the Time  1.18855E-05 7.49265E-06 2 Variance 3976.00 3975.38 3975.50 Thresholds 86 155 85 155 85 154 Time(s) 0.305679 0.145634 0.170387 Variance of the Time  1.24418E-04 6.29941E-06 3 Variance 4113.80 4112.05 4112.05 Thresholds 79 141 195 80 141 195 80 141 194 Time(s) 18.482669 2.997714 0.247683 Variance of the Time  7.63218E+01 1.6896E-05 4 Variance 4.1825e+03 4181.03 4181.56 Thresholds 67 110 158 199 67 109 158 200 66 110 157 198 Time(s) 1260.671565 0.647337 0.704557 Variance of the Time  1.27320E-01 4.11583E-02Table 4
Maximum entropy multilevel thresholding.
Image Number of thresholds Index Exhaustion method PSO SSEO Female 1 Entropy 8.8020 8.8015 8.8020 Thresholds 101 100 100 Time(s) 0.102574 0.112157 0.127748 Variance of the Time  5.25562E-05 9.34983E-07 2 Entropy 12.4804 12.4794 12.4803 Thresholds 78 154 78 153 77 153 Time(s) 1.260250 0.269593 0.227045 Variance of the Time  2.69640E-02 6.11216E-06 3 Entropy 15.6532 15.6440 15.6452 Thresholds 56 105 157 55 105 158 56 104 160 Time(s) 96.463659 0.465068 0.392787 Variance of the Time  1.87901E-02 1.20258E-05 4 Entropy 18.5869 18.5527 18.5644 Thresholds 55 103 151 189 56 103 151 190 55 103 150 188 Time(s) 8569.387408 1.410479 1.202857 Variance of the Time  1.38343E+00 2.44926E-04  Jelly beans 1 Entropy 8.0707 8.0698 8.0707 Thresholds 165 164 164 Time(s) 0.100219 0.111406 0.125397 Variance of the Time  9.62990E-05 5.44680E-06 2 Entropy 11.5128 11.5097 11.5127 Thresholds 104 167 104 166 103 166 Time(s) 1.179629 3.592932 0.214063 Variance of the Time  5.54159E+01 4.61903E-06 3 Entropy 14.4455 14.4372 14.4398 Thresholds 183 126 169 83 126 167 82 123 167 Time(s) 98.183702 0.986705 0.395346 Variance of the Time  1.40280E+00 6.48835E-06 4 Entropy 17.0383 17.0220 17.0224 Thresholds 63 94 131 169 62 94 131 167 62 93 130 167 Time(s) 5906.410023 4.553701 1.213194 Variance of the Time  3.14459E+01 2.5204E-05  San Diego 1 Entropy 8.9271 8.9270 8.9270 Thresholds 135 134 134 Time(s) 0.109754 0.112083 0.141757 Variance of the Time  1.14090E-04 2.89662E-06 2 Entropy 12.4042 12.4016 12.4041 Thresholds 60 141 59 141 59 140 Time(s) 1.227037 0.237745 0.244125 Variance of the Time  5.27638E-03 1.02371E-05 3 Entropy 15.5268 15.5216 15.5243 Thresholds 57 113 164 57 112 164 55 110 163 Time(s) 99.006925 0.453377 0.419603 Variance of the Time  7.76216E-02 7.25335E-06 4 Entropy 18.3835 18.3722 18.3729 Thresholds 54 101 143 183 54 101 142 183 53 100 142 181 Time(s) 6188.493106 2.247213 1.270564 Variance of the Time  8.36189E+00 2.08382E-05  Sailboat 1 Entropy 8.9890 8.9890 8.9890 Thresholds 125 125 124 Time(s) 0.112816 0.117312 0.140577 Variance of the Time  3.41402E-05 5.37697E-06 2 Entropy 12.5254 12.5242 12.5254 Thresholds 92 163 92 164 91 162 Time(s) 1.257669 0.210157 0.242673 Variance of the Time  3.78777E-03 7.77268E-06 3 Entropy 15.5654 15.5611 15.5619 Thresholds 73 120 170 73 120 169 73 120 170 Time(s) 101.733913 0.609797 0.420331 Variance of the Time  9.89256E-02 8.98313E-06 4 Entropy 18.3682 18.3577 18.3582 Thresholds 71 113 157 196 71 112 157 195 69 111 156 195 Time(s) 6206.751312 2.575215 1.287299 Variance of the Time  7.95833E+00 2.50911E-04From the data of the tables, in both Otsu maximum interclass variance (see Table3) and Kapur maximum entropy (see Table 4) thresholding, the exhaustion method is fastest in one threshold situation for all the four images because the computational complexity of this situation is not so high and the advantage of the heuristic methods is not exerted. The time efficiency of SSEO is lowest, for the additional commutating time introduced by subspace partition processing cannot be ignored compared with the whole calculation time.As the number of threshold increases, the computational complexity exponentially increases, the advantage of heuristic methods emerges, and subspace elimination effect on improving search efficiency is also protruded. Starting from 2 thresholds situation, SSEO is superior to PSO in time consuming. This means that it performs better than the optimization algorithm used in searching for best point of subspaces. This is from the optimization of search space but improvement of the algorithm itself. Other heuristic algorithms could also be combined with the presented idea to improve optimization efficiency. The deviation of calculation time of SSEO is much smaller than PSO which implies that SSEO is of more stable performance in the searching course.


## 5. Conclusions

Based on the observation that when searching something in a physical or logic area, people usually partition the whole potential area into some subspaces to search and eliminate the subspaces that the target unlikely exist in and to explore the rest again after a preliminary expedition, an optimization method based on subspace elimination strategy and heuristic algorithms is proposed and applied to multilevel thresholding problem. Its essential idea is to optimize the position to be searched, but the search method itself to enhance the search efficiency of optimization. Instead of reducing computational complexity, it improves search speed by improving the effectiveness of searching work. The real global best of the whole decision variable space is selected from the global bests of every subspace by elimination. Compared with using a heuristic algorithm in the whole space directly, using it in a subspace is beneficial to reduce the probability of entering the local minimum with the explored space of it being much smaller. In principle, any heuristic algorithm can be used in subspace optimization, and its work effectiveness is to be enhanced compared with that used directly in the whole space. Theoretical analysis and experimental results show that the presented method has the advantage of compressing optimization efficiency. In addition, because there is little relevance between the works in every subspace in the method, it is especially propitious to be used in parallel computing environment.Experiments also show that efficiency of the presented method is not as high as exhaustion in simpler optimization problems. This is mainly due to the fact that the time expended in eliminating operation takes a relatively large part of all the time in the simple computation situation. From theoretical analysis, it can be found that the number of subspaces increases exponentially with the dimension of the decision variable space. It may be anticipated that the efficiency of this method would decrease when used in optimal problem of high-dimensional space. The strategy of elimination process and subspace division mode can influence the performance of the presented method and may be further researched.



---

*Source: 6706590-2019-06-25.xml*