A Multilevel Image Thresholding Method Based on Subspace Elimination Optimization

the original


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][2][3][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 (t 1 , t 2 , ⋅ ⋅ ⋅ t k ) 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][10][11].Generally, heuristic algorithms are inspired by such things as natural phenomena, physical laws, biological social activities, and evolutionary processes [10][11][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][14][15].These optimization methods combined with Otsu's and Kapur's segmentation index are widely used in multithresholding segmentation [16][17][18][19][20][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 2 Mathematical Problems in Engineering 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][31][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][35][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.

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.

Subspace Division and Identification.
For a function defined in n-dimensional space f(), how to find the optimum of f() is its optimization problem.In a digital quantity (such as grey level of digital image) space, decision variable  can be defined as: x j ∈ {0, 1, 2, . . .
is a positive integer.For example, when x j express the grey level of a pixel of digital image with 8-bit binary number,   = 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 in j th dimension is marked as   .Each

Subspace label
Region in dimension 0 Region in dimension 1 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(k 0 , k 1 , ⋅ ⋅ ⋅ k −1 ), in which k j is the segment order number in jth dimension.The region of B(k Figure 1 shows that the examples of dimension n are, respectively, 1 and 2, and every dimension is partitioned into 4 segments.

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): where   ( + 1) is component of the position of particle i in dimension j at time ( + 1),   (),   () 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,  1 and  2 are coefficients,  1 () and  2 () are random number in [0, 1], and  = 1, 2, ⋅ ⋅ ⋅   ,  = 0, 2, ⋅ ⋅ ⋅  − 1.

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 number   , search times   , search round number t=0, and the state of all subspace blocks to "active." Step 2. Evenly distribute particles in every subspace •   (6) where P is jth dimension component of the position of the particle  at t.
Step 3. Calculate the initial optimum value of B(k Step 4. Determine the initial best position in B(k Step 5. Search all the "active" subspace blocks   times with PSO by   particles.
Step 6. Assign the new optimum value of B(k . Step 7. Record the coordinate of the new best position in B(k 0 , k 1 , ⋅ ⋅ ⋅ k −1 ): Step 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 Step 11.
Step 10.Adjust the subspace block search parameters   and   ; 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.

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: where 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, 255]x[0, 255]x[0, 255].If the coordinate of the optimum point is (t 1 , t 2 , t 3 ), the three thresholds are t 1 , t 2 , t 3 , respectively.
For an L grey-level digital image G which is of N  column and N  row, its pixel number is N  xN  , and grey-level set is {0, 1, 2,  − 1}.The probability with a pixel having a grey level g  is 11) where pixels in image G are divided into k + 1 classes C 0 , C 1 , . .., and C k according to their grey level with k thresholds The probability   of a pixel belonging to the class C i and the grey level probability mean   of class C i can be calculated as follows: where i = 0, 1, 2, . . ., k, t 0 = 0, and t k+1 = L.The grey level mean value of the image is The Otsu interclass variance [5] is The Kapur entropy value [6] is The Otsu maximum between class variance multilevel thresholding is to find the thresholds t 1 , t 2 ⋅ ⋅ ⋅ t  with which  2 is maximized, and the Kapur maximum entropy multilevel thresholding is to find the thresholds t 1 , t 2 ⋅ ⋅ ⋅ t  with which Ψ is maximized.

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 Figure 2. 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 Win-dows10.Software platform through which the algorithm is programed is MATLAB R2016a.
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 parameters  1 = 2,  2 = 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, 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.
From the data of the tables, in both Otsu maximum interclass variance (see Table 3) 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.

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.

Data Availability
The original images used in the experiment as illustrated in Figure 2 are available at the USC-SIPI image library.The parameters about the experiment and the results of the experiments to support the findings of this study are included within the article as Tables 2, 3, and 4. The program (in MATLAB only) of the algorithm used to support the findings of this study is available from the corresponding author upon request.

Figure 1 :
Figure 1: Examples of subspace division and identification.

Table 1 :
Subspace division and identification.