The Evaluation of Nature-Inspired Optimization Techniques for Contrast Enhancement in Images: A Novel Software Tool

This paper primarily focuses on developing an application to explore the integration of optimization theory within the multi-parametric contrast-based IE framework. In this regard, the HS technique was selected for contrast enhancement, and the optimization techniques were applied to determine the tuning parameters. In addition to a parameter-dependent technique, a non-parametric contrast enhancement technique, Automatic Histogram Equalization, was included in the study for comparative analysis of the results.

3.1 Employed dataset

This study utilized images already distorted from the Tampere Image Dataset (TID2008) [31] for tests, comparing the results with those obtained using the standard non-parametric histogram equalization technique. The TID2008 dataset, curated by Ponomarenko et al., encompasses 25 reference images subjected to 17 different distortion effects across four levels. This compilation amounts to a total of 1700 images, each with a resolution of 512×384 pixels. Given that HS, and equalization predominantly affect pixel values, the focus was on adjusting contrast levels. Therefore, only images with contrast-based distortions and their respective levels were selected. The dataset incorporated 100 distorted and 25 original images to evaluate the performance of various techniques. An example of an image subjected to different levels of contrast-based distortion is presented in Figure 1. HS with parameters optimized by different techniques was also tested on original images. This testing aimed to validate the proposed software concerning its impact on distortion. The performance was quantitatively assessed based on several metrics.

1.png

Figure 1. a) Original Image, b) Contrast Increment +1, c) Contrast decrements -1, d) Contrast increment +2, e) Contrast decrements -2

3.2 Standard histogram equalization

Histogram equalization represents one of the IE techniques used for adjusting contrast levels. This method involves the automatic alteration of the histogram curve using the pixel-based color space present in the image [32]. The process begins with the computation of the probability mass function (PMF) for all pixels within the image. Subsequently, the cumulative distribution function (CDF) is calculated utilizing the PMF and multiplied by various levels to determine new pixel intensities. The respective formula is presented in Eq. (1), where L symbolizes the multiplication level for the computed CDF and p_n denotes the bit number of possible color or intensity values. The terms f_i.j and g_i.j represent the coordinates of pixels corresponding to the i^th row and the j^th column, and the generated new image pixel, respectively.

$g_{i . j}=f l o o r\left((L-1) \sum_{n=0}^{f_{i . j}} p_n\right)$                 (1)

3.3 Histogram stretching

HS is a simple and effective technique for improving the contrast in an image. It works by spreading out the most frequent intensity values or colors in an image, thus enhancing the contrast. The key parameters of HS are the minimum and maximum intensity values of the image. By stretching the range of intensity values, the contrast of the image can be increased. In the HS technique different mathematical operations can be utilized to enhance the contrast levels [33]. We selected multiplication, addition, and root extraction mathematical operations to be used in HS technique.

The parameters within these operations serve as vital components for obtaining an efficient enhancement algorithm. It is crucial to understand that these parameters are image-specific, meaning a singular parameter cannot be universally applied to all images. Therefore, an adaptive process for parameter adjustment is necessary for each distinct image. In response to this requirement, a software tool incorporating nine nature-inspired optimization techniques for optimizing defined HS parameters is presented within this study. The results procured by this software are thoroughly discussed.

2.png

Figure 2. Utilized nature-inspired optimization techniques

The techniques employed within this study are itemized in Figure 2. Each optimization process necessitates a convergence point to terminate the iteration and define the optimal parameters. For this study, Eq. (2) is used to establish these convergence points.

$S_t=\left(\operatorname{argmax}\left(I_t\right)-\operatorname{argmin}\left(I_t\right)+1\right) \times n_t$              (2)

where, I represents the intensity values of image t, and n refers to the number of colors in the histogram chart of image t. The invariability of S within a certain range serves as the convergence criterion for the employed optimization techniques. If the S value remains stable during the last k iterations of the optimization process, then the HS with optimization techniques will be terminated. k can be selected as an arbitrary constant. The stopping criterion is given in Eq. (3).

$\begin{aligned} S_t-20 & <S_t<S_t+20 \\ S_t & =\frac{1}{k} \sum_{q=1}^k S_t^q\end{aligned}$              (3)

Table 1. The key properties of each optimization technique, including the type of problems it is best suited for, its complexity, the number of parameters it requires, its advantages, and its limitations

$g_{i . j}=\sqrt[\gamma]{\alpha \times I_{i . j}+\beta}$             (4)

where, I_i.j is the intensity values of image at i^th and j^th location. α, β and γ serve as the tuning parameters of the stretching as scaling, shifting, and degrading of the image contrast, respectively. Optimization methods will be applied via the proposed software for the selection process of α, β and γ parameters. g_i.j indicates the new pixels obtained through HS. The nature-inspired optimization techniques utilized in the proposed software will be briefly elucidated in the following subsections. Additionally, Table 1 summarizes the key properties of each optimization technique, including the type of problems it is best suited for, its complexity, the number of parameters it requires, its advantages, and its limitations. This will provide a clear overview and facilitate a straightforward comparison of the different optimization techniques.

3.3.1 Particle swarm optimization (PSO)

In the PSO algorithm, a "swarm" is composed of combinations of particles, where each element represents a "particle". The particles initiate the exploration for the optimal solution by taking random values in the solution space. Each particle is composed of two vector components: (a) the position vector 'x', and (b) the velocity vector 'v'. The position vector contains the positional information of the particle, while the velocity vector maintains the information about the change in position and direction of the particles. The essence of PSO lies in transitioning each particle's current position in the swarm to the best personal position ('pbest') identified previously, and the best global position ('gbest') discovered so far. All particles aim to orient themselves towards these two best positions using their velocity vectors. The updating process of velocity and position is calculated using Eq. (5).

$\begin{gathered}x_i^k=(\alpha \beta \gamma) \\ v_i^{k+1}=W . v_i^k+c_1 . \operatorname{rand}_1^k .\left(\text { pbest }_i^k-x_i^k\right)+c_2 . \text { rand }_2^k .\left(\text { gbest }^k-x_i^k\right) \\ x_i^{k+1}=x_i^k+v_i^{k+1}\end{gathered}$              (5)

where, c₁ and c₂ are the constants for learning parameters leading the particles to its pbest and gbest position. In this study, 2 is selected for $c_1$ and $c_2$. $\operatorname{rand}_1^k$ and $\operatorname{rand}_2^k$ are randomly selected numbers between 0 and 1 in uniform distribution at iteration k. W represents the weight of inertia to be used in definition of balance between local and global searching. v_i^k and x^k_i indicates the velocity and position of i_th particle in swarm at iteration k. x refers to the vector including α, β and γ parameters which needs to be optimized in the HS equation Eq. (4) [9-11, 34]. PSO optimizes the HS parameters by mimicking the social behavior of bird flocking. Each particle represents a potential solution and adjusts its position in the search space based on its own best experience and the best experience of the entire swarm. While PSO can quickly converge to an optimal solution, it is sensitive to initial conditions and can suffer from premature convergence.

3.3.2 Artificial bee colony (ABC)

The ABC algorithm comprises two primary components: bees and a food source. Mirroring natural behavior, bees in this algorithm seek out rich food sources in proximity to their hive. Essentially, the ABC algorithm leverages this phenomenon. The algorithm includes three groups of bees: worker bees assigned to specific food sources, onlooker bees observing the dance of worker bees within the hive to select a food source, and scout bees randomly searching for food sources. The ABC algorithm is initialized with Eq. (6) for the HS problem [12-14, 35].

$x_i=l_i+\operatorname{rand}(0.1) .\left(u_i-l_i\right)$             (6)

where, u and l indicate the upper and lower boundaries of the dimension i. rand is a parameter similar to that in the PSO algorithm, representing a randomly selected number between 0 and 1 with a uniform distribution. x stands for the i^thsolution in the swarm. After first iteration, each worker bee (x_i) generates a candidate solution (v_i) by using Eq. (7).

$v_i=x_i+\Phi i \times\left(x_i-x_j\right)$             (7)

x_j is the randomly selected worker bee, under the condition that i and j must not be the same. φ_i is another randomly selected number within [-1, +1] to weight the current worker bee to the optimal solution. If the candidate solution which is represented as v_i is better than its parent x_i, then update current bee to the candidate by φ_i. The best values of v_i used as α, β and γ parameters in the HS formula as given in Eq. (8). In terms of HS, each food source represents a possible solution, and the nectar amount corresponds to the quality (fitness) of that solution. While ABC has fewer control parameters and is simple to implement, it can be slow to converge in complex problems.

$v_i=[\alpha \beta \gamma]$                (8)

3.3.3 Genetic Algorithm (GA)

GA adopts the principle of evolution, favoring the survival of the fittest individuals in a population and removing the weaker ones. In GA, each individual within the population is represented as a chromosome. The total quantity of chromosomes indicates the size of the population. Each chromosome carries specific values, reflecting their fitness according to the employed GA fitness function. Initially, these chromosome values are selected at random. Subsequently, these values either increase or decrease in accordance with the utilized fitness function. Values that increase gain importance for the next iteration as they come into play during the crossover process. The crossover mechanism allows for superior individuals to intermingle and generate better offspring within the population. Following this, chromosomes undergo mutation at a predetermined rate to avoid becoming stuck in local minima or maxima in the problem space. The algorithm repeats these steps until a predefined stopping condition is reached. The final chromosome, discovered at the conclusion of all iterations, yields the optimal solution to the problem. In this setup, chromosomes comprise three values: α, β, and γ. The optimal values are applied to the HS of an image in accordance with a specific formula Eq. (4). The fitness value is defined as per Eq. (2), following the roulette wheel selection method [25, 26]. In this framework, the HS parameters are treated as genes in a chromosome, and a population of these chromosomes evolves over generations through crossover, mutation, and selection operations. Despite offering global search capabilities, GA can be computationally demanding and necessitate careful parameter tuning.

3.3.4 Differential evolution (DE)

DE algorithm is a simple but powerful algorithm that works population-based like GA. During the iterations in the algorithm, better results are searched for the solution of the problem with the help of operators- selection, mutation, and recombination. In contrast to binary GA, DE uses variables with real values so it can be implemented on continuous problems. Additionally, every operator is not sequentially applied to the entire population. In terms of mutation operator, the difference between GA and DE can be explained as the result of small perturbations to the genes of an individual and the result of arithmetic combinations of individuals, respectively. Especially, the mutation operator is adjusted progressively and updated evolutionary to achieve the best value. In other words, it is not specified as predefined function [15-18, 36]. In HS, DE can identify the optimal set of stretching parameters by balancing exploration and exploitation. However, DE can be sensitive to parameter settings and may suffer from slow convergence speed in high dimensional problems.

3.3.5 Simulated Annealing (SA)

SA is designed to find the maximum or minimum values of functions with a large number of variables. This algorithm and its derivatives have been inspired by the similarity between searching for solutions to an optimization problem and the process of annealing in metallurgy. Annealing describes the process of heating a metal to a certain degree and then gradually cooling it. The annealing process starts with a high temperature value. At this temperature, a random solution is selected from the solution space and the neighbors of this solution are examined by the number of defined iterations. If any of the neighboring solutions produces a better fitness value, then the process is continued with this neighboring solution. In case of failure in neighboring solutions, value is assigned randomly in order to avoid from the local maximum and minimum. When the temperature decreases, the probability of random selection of solutions will also decrease. This probability is calculated according to the Eq. (9).

$\begin{aligned} T & =T \times a \\ \mathrm{~N} & =\mathrm{N} \times \mathrm{k}\end{aligned}$                (9)

The temperature and the number of iterations are reduced according to Eq. (9) where a and k are the predetermined constants in the range (0, 1). T is the temperature and N is the number of iterations. When the temperature reaches its minimum value, the algorithm terminates and determines the optimal α, β and γ parameters [19, 20]. In HS, SA can explore the solution space thoroughly. But it requires careful temperature scheduling and can be computationally intensive.

3.3.6 Bacterial foraging optimization (BFO)

Bacterial Foraging Optimization (BFO) offers a model for the food search behavior of bacteria and simulates certain movements such as "tumbling" or "swimming". BFO's operation can be examined in four processes: Chemotaxis, Swarming, Reproduction, and Elimination. In Chemotaxis, two movements are defined: tumbling and swimming. Swimming refers to moving in the same direction for a certain number of steps, while tumbling represents the random movements of bacteria in search of food sources. In Swarming, bacteria move collectively, forming concentric groups towards the food sources. During reproduction, unfit bacteria are eliminated and those that are sufficiently fit reproduce by splitting. In the elimination step, several bacteria can be eradicated or relocated due to sudden environmental changes. These steps are mathematically implemented in the BFO optimization technique, which was employed in the proposed software to determine the optimal α, β, and γ values. For HS, BFO is capable of thoroughly searching the solution space. However, it should be noted that BFO can suffer from slow convergence speed and high computational costs [21, 22, 37].

3.3.7 Cat swarm optimization (CSO)

CSO is inspired by the behavior of cats, particularly their food-seeking behavior. CSO incorporates two main steps: Seeking and Tracing. In seeking mode, parameters represent the state of a cat visually scanning its environment to decide on its next destination. In tracing mode, parameters depict a cat as it tracks a target and moves closer to it [23, 24, 38, 39]. In the context of HS, it can adaptively select the appropriate mode to locate the optimal stretching parameters. However, its performance may be sensitive to the proportion of cats in each mode and the number of cats, which are parameters that need to be carefully tuned.

3.3.8 Dragonfly Algorithm (DA)

The DA draws its main inspiration from the static and dynamic swarm behavior of dragonflies in nature. The two primary stages of optimization-exploration and exploitation-are modelled after the social interactions of dragonflies during navigation, food search, and evading enemies in static and dynamic scenarios. During the exploration phase, dragonflies form sub-swarms within static swarms, and each of these sub-swarms explores different areas. In the exploitation phase, these sub-swarms communicate with each other and converge to move in one direction within the larger static swarm. From a programming perspective, the position information of dragonflies begins with a certain number of randomly generated particles. The satisfaction rates for all of these particles are subsequently calculated. Distances to food and enemies are calculated based on these satisfaction rates. The parameters are updated according to the dragonflies’ movements, positions, fitness values, and distances to food and enemies using Eq. (10) [27, 28, 40].

$\begin{gathered}S_i=\sum_{j=1}^N X-X_j \\ A_i=\frac{\sum_{j=1}^N V_j}{N} \\ \mathrm{C}_{\mathrm{i}}=\frac{\sum_{j=1}^N X_j}{\mathrm{~N}}-X \\ F_i=X^{+}-X \\ E_i=X^{-}+X\end{gathered}$              (10)

where, S refers to distribution of dragonflies and equals to the sum of the distances between the dragonfly and its neighbors. A indicates the harmony of dragonflies and is equal to the average speed of flocks. C is the distance of each dragonfly to the average positions of flocks and expresses the dragonfly conjunction. F states the distance between the dragonfly and the food source. E refers to the enemy’s distraction and equals to the sum of the position of the dragonfly and the enemy. Lastly, W indicates immobility. These parameters are updated at each step, and the position vectors of the dragonflies, which refers to α, β and γ in our study, are updated by using Eq. (11). For HS, DA provides an efficient exploration and exploitation mechanism. However, parameter tuning is crucial in DA and it can get stuck in local optima.

$\begin{gathered}\Delta \mathrm{X}_i^{\mathrm{t}+1}=\left(\mathrm{sS}_i+\mathrm{aA}_i+\mathrm{cC}_i+\mathrm{fF}_i+\mathrm{eE}_i\right)+\mathrm{w} \Delta \mathrm{X}_i^{\mathrm{t}} \\ X_i^{t+1}=X_i^t+\Delta X_i^{t+1}\end{gathered}$                (11)

3.3.9 Black hole algorithm (BH)

In the BH, the best candidate is chosen as the “black hole” among all candidates in each iteration by using a fitness function. All other candidates are assigned as regular “stars”. Then stars are drifted to the black hole based on their current location and attraction force. If a star gets too close to the black hole, it is swallowed by the black hole and disappears. In this case, a new star is randomly generated, placed in the search field, and a new search is started. The implemented fitness and update functions are given in Eq. (12) and Eq. (13), respectively. Similar to other optimization techniques, updated state parameter (X_i(t)) indicates the optimized parameters α, β and γ in this study [29, 30]. BH excels in global optimization problems but might struggle with local optimization, as it primarily focuses on the global best solution. The performance of BH can be greatly influenced by its only parameter, the event horizon.

$f_{B H}=\sum_{i=1}^{\text {pop_size }} \operatorname{eval}(p(t))$                  (12)

$X_i(t)=X_i(t)+\operatorname{rand} \times\left(X_{B H}-X_i(t)\right)$                (13)

3.4 Performance evaluation

Enhanced versions of the distorted images were compared with reference images using various optimization techniques for HS, as well as standard histogram equalization techniques. Performance evaluation was carried out using certain image quality metrics such as Mean Square Error (MSE), Peak Signal to Noise Ratio (PSNR), and Structural Similarity Index (SSIM) [41]. MSE was computed as described in Eq. (14).

$M S E=\frac{1}{N \times M} \sum_{i=0}^{N-1} \sum_{j=0}^{M-1}[X(i . j)-Y(i . j)]^2$              (14)

where, N and M represent the total number of pixels in images as width and height size, respectively. X (i.j) and Y (i.j) represents the pixel values at i_th row and j_th column of original and contrast enhanced of the distorted images. The PSNR metric, derived from the MSE, can be calculated using Eq. (15).

$P S N R=20 \log _{10}\left(\frac{255}{\sqrt{M S E}}\right)$               (15)

A higher PSNR value indicates greater symmetry in images according to the PSNR metric. In the case of identical images, PSNR is infinite or undefined due to the MSE value being "0". SSIM is another image quality measurement that primarily quantifies image quality degradation. SSIM constructs perception-based models that regard image degradation as perceived change in structural information. The formula for SSIM is given in Eq. (16).

$\operatorname{SSIM}(x . y)=\frac{\left(2 \mu_x \mu_y+c_1\right) \times\left(2 \sigma_{x y}+c_2\right)}{\left(\mu_x^2+\mu_y^2+c_1\right) \times\left(\sigma_x^2+\sigma_y^2+c_2\right)}$            (16)

where, µ_x, µ_y, σ_x, σ_y and σ_xy are the local means, standard deviations, and cross-covariance for images x and y, respectively.

Figure 3. The flowchart of the experiments