{'fraction_non_alphanumeric': 0.0379124175164967, 'fraction_numerical': 0.015436912617476505, 'mean_word_length': 4.511571523032841, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multilevel Thresholding of Brain Tumor MRI Images: Patch-Levy Bees Algorithm versus Harmony Search Algorithm

Image segmentation of brain magnetic resonance imaging (MRI) plays a crucial role among radiologists in terms of diagnosing brain disease. Parts of the brain such as white matter, gray matter and cerebrospinal fluids (CFS), have to be clearly determined by the radiologist during the process of brain abnormalities detection. Manual segmentation is grueling and may be prone to error, which can in turn affect the result of the diagnosis. Nature inspired metaheuristic algorithms such as Harmony Search (HS), which was successfully applied in multilevel thresholding for brain tumor segmentation instead of the Patch-Levy Bees algorithm (PLBA). Even though the PLBA is one powerful multilevel thresholding, it has not been applied to brain tumor segmentation. This paper focuses on a comparative study of the PLBA and HS for brain tumor segmentation. The test dataset consisting of nine images was collected from the Tuanku Muhriz UKM Hospital (HCTM). As for the result, it shows that the PLBA has significantly outperformed HS. The performance of both algorithms is evaluated in terms of solution quality and stability.


INTRODUCTION
Cancer is the second mortality causing factor in the world with 8.8 million cases reported as of 2015 [1]. In Malaysia, cancer is the third highest contributor to mortality [2]. Categorically speaking, Malays have the highest survival rate among those affected by brain and nervous system cancers, which amounts to 43%, as reported by the Malaysian Study on Cancer Survival (MySCan), 2018. This has become the second highest rate in terms of survivors, with ovarian cancer (54.8%) third, respectively [1]. Indians and Chinese have the highest survival percentage in other types of cancer [1]. It can be concluded that Malays have a higher risk of developing brain cancer in Malaysia.
Brain tumor or brain cancer is a cluster of abnormal cells in the human brain [3]. There are two types of tumors, namely malignant (cancerous) and benign (non-cancerous) [3]. Medical imaging is an important medium for radiologists to accurately diagnose brain disease such as cancer. High resolution MRI of the brain is needed for better brain tumor detection.
The advantage of MRI is that it is the least risky method for constructing data with spatial resolution in high scale and non-invasive modality compared to other techniques of diagnostic imaging [4]. There are several methods that have been applied successfully in MR image segmentation.
Manual segmentation of MRI images is arduous, expensive and time consuming. Any mistake is susceptible due to blurriness on tissues boundaries, low tissue contrast and poor hand-eye coordination. Consequently, discrepancies are common among radiologists determining a variety of structures. Therefore, a plethora research studies were introduced to increase the accuracy of brain tissues segmentation upon brain disease diagnosis. Effective brain MRI segmentation can improve a classification of brain diseases with better precision. Joliot et al. [5] proposed a technique of three-dimensional automatic segmentation that separates gray level thresholding of gray matter and white matter. In fact, brain MRI images occasionally have some noise component as a result of patient movement during brain scanning. Raquel et al. [6] introduced a hybrid model for multispectral brain MRI segmentation based on a radial basic network. This method falls short of solving the brain MRI images with noise. To overcome this problem, a clustering method is put forth.
In a number of pixel-based approaches, the clustering method has the most prominent application in terms of image segmentation [7]. In the last decade, clustering-based approaches garnered significant interest in the domain of medical imaging. The fuzzy cmeans (FCM) clustering algorithm was introduced by Bezdek [8] to minimize an objective function which is known as fuzzy membership and a set of cluster centroids. However, FCM still blunders in separating the pixel since it contains noisy pixels.
Histogram-based thresholding is a very popular tool in image segmentation. Bi-level thresholding is also known as the simplest problem compared to multi-level thresholding. Bi-level thresholding segments the image into two groups. If the image consisting of a gray-pixel level is higher than the threshold, it is defined as one group. Otherwise, it is included as another group. In the case of multilevel thresholding, it is more challenging to define a group of pixels when more details of segmentation are produced. Determining distinct valleys in a multilevel histogram is no easy task. Hence, the problem of multilevel thresholding has become the point of interest among researchers.
Otsu's method [9], which is also known as a nonparametric approach, selects optimal thresholds by maximizing the between-class variance of gray levels. Gray levels of the image are normally distributed [10]. This method is easy and powerful in bi-level thresholding [10]. Otsu's method can possibly be applied in multilevel thresholding. However, it is significantly formidable to determine optimal thresholds due to exponential growth in computation time. Owing to that, numerous methods for solving the multilevel thresholding problem have been proposed [11].
Multilevel thresholding-based metaheuristics is proposed by the researchers to increase search speed as it has been attested to yield the best result in (optimal) thresholds. Brain MRI segmentation focuses on the normal brain by using bacterial foraging optimization (BFO) [12], a genetic algorithm (GA) [13], and adaptive bacterial foraging optimization (ABF) [14], which have been successfully implemented.
Brain tumor MRI images encompass various artifacts such as cyst, solid and edema. Multilevel thresholding of brain tumor MRI images segmentation is also successfully applied based on metaheuristic methods such as harmony search (HS) [15]. The advantage of the Patch-Levy Bees Algorithm (PLBA) is that it is more realistic and it is attained by speeding up the search process to find an optimal threshold based on multilevel thresholding for image segmentation [16]. The PLBA has never been applied in brain MRI segmentation. Owing to that, compared to HS, the PLBA is deemed interesting for the application in brain tumor MRI images, by using the recent methods to provide better brain tumor image segmentation.

IMAGE THRESHOLDING
Image thresholding is one of the types of image segmentation techniques. It segments an image into different intensity levels of color or gray images. Image segmentation through image thresholding is the most convenient and effortless method effortless. Threshold values are the intensity values selected from the color or gray intensities of the image. Based on selection of these threshold values, an image is segmented into different regions.

GLOBAL THRESHOLdING
Global thresholding has one threshold to segment into two partitions [17]. It operates in such a way that if any gray value of images is higher than the threshold value, it will join into one group; otherwise, it will join into the other group.
Even though global thresholding is so easy to implement, it only considers segmentation of the image into two partitions. It may suffer from the loss of sharp pixels to classify in the image whose color intensities are more similar to other uniquely classified pixels. The global thresholding method is also known as a bi-level thresholding method.

MuLTILEvEL THRESHOLdING
Multilevel thresholding considers segmentation of more than one threshold value into several groups [17]. The number of segmented images is one more than the number of threshold values chosen. It can be concluded that if three threshold values are selected, it will produce four groups of regions for the segmented image. Generally, some of the sharp pixels that went through global thresholding may be lost or missing after performing segmentation separately. For the purpose of solving the shortcomings of global thresholding, multilevel thresholds are selected.
Gray images have L levels of gray in the range (0, 1, 2, ... , L-1 1). Let h(1), h(2),..., h(L-1), be the observed graylevel frequency. The probability of occurrence of the gray level i can be defined as in (1): where h(i) is the number of pixels corresponding to the gray level i and N is the total number of pixel in the image, which is equal to

HARMONy SEARCH
Harmony search (HS) was introduced by Geem in 2001 [18], based upon the idea that it can imitate the behavior that a musician needs, a well-balanced harmony. The advantage of HS is that it has a simple model and good performance in the global search domain. Furthermore, HS has been applied in numerous real-world optimization problems [19][20][21]. HS steps are shown in Algorithm 1.
The first step of HS is to initialize the harmony memory size HMS, harmony memory consideration rate HMCR, maximum dimension D, the pitch adjusting rate PAR and the bandwidth BW.
Then, the second step is to initialize harmony memory x, as in (2): (2) where x j,L and x j,U are lower bounds for a decision variable. In step 3, check the memory consideration of HS. If R 1 is smaller than HMCR, execute as in (3): where x r,j is the index random of the harmony memory for one dimension, r∈1,2,3,…HMS and j∈{1,2,3,…D}, while D is the maximum number of dimensions.
After that, check the pitch adjustment condition; if R 2 is smaller than PAR and R 3 is smaller than 0.5, execute as in (4): where φ is a random number uniformly distributed in the range [-1,1] and BW is a constant value. If R 3 is larger than 0.5, execute as in (5): If R 1 is larger than HMCR, follow (6): where rand is a random number uniformly distributed in the range [0,1].
In step 4, if f(x new ), which is an objective function of the new harmony memory x new , is smaller than f(x worst ), which is an objective function of the worst harmony memory x worst , execute as in (7): Algorithm 1: HS Input: Harmony memory size HMS, harmony memory consideration rate, HMCR maximum dimension D, maximum generation G, pitch adjusting rate, PAR and bandwidth, BW Output: worst harmony, x worst , function objective of worst harmony f(x worst ) Step 1: Initialize HMS, HMCR, PAR, D, G, and PAR.
Step 2: Initialize harmony memory HM as below: For j=1 to D Execute as in (2).
End For Calculate f(x i ) End For Step 3: Improvise g=0, while the stopping criterion is inadequate (g< G) do For j=1 to D do If R 1 ≤ HMCR % memory consideration Execute as in (3).
Execute as in (4). Else Execute as in (5). End If End If Else % random selection Execute as in (6).

End If End For
Step 4: Update the worst harmony.
Select the worst harmony vector x worst in the current harmony memory and calculate f( Execute as in (7). End If g=g+1 ; End While Step 5: Algorithm stops when the best solution is obtained.

PATCH-LEvy BEES ALGORITHM (PLBA)
The PLBA algorithm was proposed by Wasim et al. in 2016 [16]. During initialization, local search and global search of the PLBA, a patch concept and the Levy distribution were modeled. Using patch with the PLBA helps us distribute solutions in the search space, while the Levy distribution is used to get the benefit of its concept and step size. Levy distribution step size usually comprises several short steps followed by occasional longer jumps. Short steps make the exploitation of the most promising area on which the majority of the population is focused. However, a rarely long jump keeps exploring the search space by part of the population. Generally, patch and levy flights make the PLBA track more than one pick in terms of the maximization problem. Using a suitable step size of levy flights in initialization and global search, the PLBA maintained the diversity of the solution because of the rarely long jump.
The PLBA is attested to select accurate thresholds for image segmentation [16]. In comparison to Kapur's objective functions, Otsu's show a better segmentation result.
The PLBA solution represents a combination of r thresholds. Hence, the bee solution can be represented as in (8). (8) where t 1 ∈[0,L 1 ] and t i <t i+1 for all t∈ [1,r].
The PLBA-based algorithm steps are shown in Algorithm 2. Generally, initializing the threshold values is done with levy step γ=γ 1 as in (9), while in global search the levy step γ=γ 3 . Local search is performed as in (10).
where c j is the center of the j th patch, r ∈ uniform(0,1), and gives the direction to the flight.

OTSu'S CRITERION
Otsu's thresholding is a well-known method in image segmentation. It is used to choose optimal thresholding by maximizing the variance between different classes and minimizing it within the same class.
Let h be the gray-level histogram of an image with L gray levels and N pixels. So h(i) is the number of gray levels i, i ∈[0,L] in the image. Then, the probability of gray level i in the image can be calculated as in (11). Step 1: Initialize each scout bees as in (8); with a combination of threshold values defined as in (9) Step 2: Evaluate the fitness of each combination of thresholds as in (18) for the PLBA-Otsu algorithm.
Step 3: Repeat for s or until a given stopping criterion is met 3.1: Select the best m; where m<n; combinations of thresholds from the population 3.2: Select the best e ; where e<m; combinations of thresholds from the population 3.3: Start local search around the best m combinations by generating more combinations of thresholds (nep) around e combinations and (nsp) threshold combinations generated around the remaining m-e combinations. Each combination of thresholds (t i ) is represented as in (8) and generated as in (10), where t bestcur is the current best combination of thresholds in the neighborhood,

3.4:
Start global search by redistributing the remaining n-m bees by creating a new solution represented as in (8) and (9) 3.5: Evaluate the fitness of each combination of the n-m combinations as in (18).

3.6: Choose the best threshold combination of local and global search
Step 4: A combination of thresholds with the best fitness value is chosen to represent the optimal thresholds.
Otsu [9] introduced three discriminant criteria based on discriminate analysis to measure separation as in (12). (12) where σ 2 B is the variance between classes, σ 2 W is the variance within class and σ 2 T is the total variance of gray levels.
In bi-level thresholding, Otsu divided the image into two classes, i.e. the first class C 1 with pixels in the range [0,t-1], and the second class C 2 with pixels in the range [t,L-1] by threshold t. The between class variance is calculated as in (13).
where w(0,t)=∑ t-1 i=0 p i is the weight of the first class C 1 ,w(t,L)=∑ L-1 i=t p i is the weight of the second class C 2 and w(0,t)+ w(t,L)=1.
The mean of the first and second classes is given as in (14) and (15).
The total mean level of the original image is defined as in (16).
Otsu's thresholding method selects the optimal thresholds t * by maximizing the between class variance of (12) as follows (17): (17) Assuming r thresholds were supposed to be chosen, two dummy thresholds t 0 =0 and t r +1=0 are used for computational convenience, and the objective function becomes as in (18):

TEST IMAGE
The MRI data for this research collected from the Hospital Counselor Tunku Muhriz (HCTM) with a sample T2 weighted of 9 axial plane brain images for 2 patients represented as Patient A and Patient B are shown in Table 1. Patient A is a 40-year old Malay male and Patient B is a 48-year old Malay female. T2-weighted images in JPEG format show water brighter and fat darker. The details of each image are shown in Table 1 and the images with their histograms are shown in Fig. 1 Uniformity measured was employed to compare optimization technique performance [22]. The misclassification error is calculated by (19): where c is the number of thresholds used to segment the image. R j is the jth segmented region and f i is the intensity level of pixels for the segmented region. μ j is the mean of the jth segmented image and N is a total of pixels in the image. Maximum and minimum intensity of an image are represented as f max -f min , respectively. Generally, misclassification errors may lie between 0 and 1.
The performance of the algorithm is better when the misclassification error is high.
Thus, the uniformity is an indication of better image quality when U=1, and of worse image quality when U=0.

B. Stability
If the performance in each run is not unique, it is evaluated by more than one run. Each run obtained different initial values. The robustness of the algorithm is based on the acceptance of its outcome. Therefore, the algorithm was run for 50 times.
The stability of the algorithm evaluated was based on the average of the mean and standard deviation for all runs. If each run hardly varies from one to another, the algorithm was considered to be stable. The mean and standard deviation are demonstrated as in (20) and (21).
where μ j is the objective function value or the fitness value at the jth run and N is the number of run.

C. Structural Similarity Index (SSIM)
SSIM measures the similarity between the original image and the thresholded image, and it is expressed as follows: where μ I and μ I~ are the mean values of the original image I and the reconstructed image I~, and σ I~ is the standard deviation of the original image I and the reconstructed image I~, σ I I~ is the cross-correlation and C1 and C2 are constants equal to 0.065. The SSIM ranges between -and +1 and the SSIM value is equal to one. This shows that the original image and the reconstructed image (a threshold image) are similar. It is presumed that the algorithm can produce acceptable p performance if the SSIM value is close to 1.

d. Peak Signal-to-Noise Ratio (PSNR)
The PSNR, which is quantified in decibels (dB), shows a dissimilarity between an input image and a threshold image. A higher value of the PSNR requires better quality of the threshold image. RMSE and the PSNR values are calculated as follows: where M × N is the size of an image, while I and J represent the pixel values of the original and decompressed images, respectively. In this experiment, researchers have taken M × N to be a square image. f(I,j) is an original image, whereas f~(I,j) is a reconstructed image.

PARAMETER SETTINGS
For PLBA-Otsu and HS-Otsu tests presented in this paper, the default parameters of all methods are listed in tables 2 and 3, respectively. The PLBA-Otsu parameter values were suggested by Wasim et al. in 2016 [16]. The Levy step size in local search (Υ 1 ) was set to a small value to support the exploitation capability of good regions, while the Levy step size in global search and initialization (Υ 2 and Υ 3 ) was set to large values to maintain the diversity of the population [16]. Basically, HS parameters are the harmony memory size (HMS), that is, the number of solution vectors lying on the harmony memory (HM), the harmony memory consideration rate (HMCR), the pitch adjusting rate (PAR), the distance bandwidth (BW), and the number of improvisations (NI), which represents the total number of iterations. These HS parameters, with the exception of NI, were obtained from Oliva et al. [20]. HS is stopped when the NI maximum is reached. Step size of Levy flights (local), (Υ 1_init ) 0.01 Step size of Levy flights (initial), (Υ 2 ) 0.5 Step size of Levy flights (global), (Υ 3 ) 0.5

ExPERIMENTAL ENvIRONMENT
Both PLBA and HS based multilevel thresholding were implemented in the Matlab Language under operation system Microsoft 10 Pro. HS used the platform of Intel (R) Core™ i7-2620M CPU @ 2.70 GHz 2.70 GHz and 4.00 GB RAM. On the other hand, the PLBA used the platform of Intel (R) Core™ i7-2640M CPU @ 2.80 GHz 2.80 GHz and 8.00 GB RAM.

SOLuTION QuALITy
The mean values of Otsu's objective functions and PLBA and HS optimal thresholds are shown in Table 4 and Table 5, respectively. PLBA and HS algorithms yield different solution quality (thresholds). In Table 4, it can be observed that the PLBA algorithms provide higher objective values than HS for Otsu's criteria in relation to all numbers of thresholds. It is triggered by the balance of PLBA search through its exploitation and exploration.
An example of segmented images and corresponding thresholds on the histogram obtained by means of PLBA and HS algorithms at threshold levels 2, 3, 4 and 5 are shown in Fig 2. In these figures, it can be observed that a segmented image has high performance when (C=5), compared with (C=4), (C=3) and (C=2).
The one-way ANOVA test is the most suitable statistical test for comparing two algorithms during performance evaluation.  The hypothesis made for the analysis of Table 6 is as follows: Null hypothesis, H0: There is no significant difference in the objective value between the two methods (the PLBA and HS).

Alternative hypothesis, H1:
There is a significant difference in the objective value between the two methods (the PLBA and HS).

Test images c
Optimal threshold values HS PLBA    The p values are greater than 0.05 (a significance level of 5%) and the F value is less than F-critical, therefore the null hypothesis cannot be rejected. Hence, it is concluded that there is no significant difference in the objective value between the PLBA and HS.

uNIFORMITy
A uniformity measure of the PLBA yields the highest values of HS (refer to Table 7) in all of the test images and thresholds, with the exception of two thresholds in slices 165, 166, 167, 168 and 38. The PLBA algorithm has successfully demonstrated and yielded the best quality of segmented images. Table 8 portrays the single-factor ANOVA test results regarding the uniformity obtained from Table 7.
The hypothesis made for the analysis of Table 8 is as follows: Null hypothesis, H0: There is no significant difference in the uniformity between the two methods, i.e. the PLBA and HS.
Alternative hypothesis, H1: There is a significant difference in the uniformity between the two methods, i.e. the PLBA and HS.
The p values are greater than 0.05 (a significance level of 5%) and the F value is less than F-critical; therefore, the null hypothesis will not be rejected. Hence, this concludes that there is no significant difference in the uniformity between the PLBA and HS.

STABILITy
Based on Table 10, the PLBA algorithm shows a lower standard deviation value for all image thresholds. Thus, the PLBA is more stable compared to HS. Table 9 presents the one-way ANOVA test results regarding uniformity obtained from Table 10.
The hypothesis made for the analysis of Table 9 is as follows: Null hypothesis, H0: There is no significant difference in the standard deviation between the two methods, i.e. the PLBA and HS.
Alternative hypothesis, H1: There is a significant difference in the standard deviation between the two methods, i.e. the PLBA and HS.
The p values are less than 0.05 (a significance level of 5%) and the F value is greater than F-critical, by virtue of which the null hypothesis is rejected. Hence, it is concluded that there is a significant difference in the standard deviation between the PLBA and HS.  Table 9. One-way ANOVA test of the PLBA and HS based on standard deviation Table 11 shows that the SSIM of the PLBA is mostly higher than HS. Some SSIM of HS is higher than the PLBA, such as slices 166, 168 and 38 for c=3, and slices 167 and 169 for c=2.

PEAK SIGNAL-TO-NOISE RATIO (PSNR)
The results of the PSNR value in all images (refer to Table 11) are directly proportional to the increasing threshold values even though HS has achieved a higher PSNR value for most images and thresholds. The PSNR value of the PLBA is the best value in c=5 for slices 9, 162, 164 and 165. Thus, well segmented images provide a better result when the level of thresholds is high.   Table 11. SSIM and PSNR

CONCLUSION
In this paper, the PLBA algorithm based multilevel thresholding using Otsu's criterion is discussed. T2 weighted brain MR axial plane images are applied as test images. It is observed that the PLBA is higher in terms of the fitness value (Otsu) and low standard deviation for all test images. The uniformity measure of PLBA based multilevel thresholding is mostly higher than HS and it represents the quality of the segmentation algorithm. The segmented image yields more information of the image when the number of threshold levels increases. The statistical performance of the PLBA and HS is evaluated by using one-way ANO-VA. The test showed that there is a significant difference in standard deviation for both algorithms even though there exists no significant difference of the objective value and uniformity between the PLBA and HS. This proves that the PLBA is more stable than HS. A higher threshold level implies the effectiveness of brain MRI images to differentiate between normal and abnormal tissues.