Infrared Image Enhancement Using Adaptive Histogram Partition and Brightness Correction

Abstract

Infrared image enhancement is a crucial pre-processing technique in intelligent urban surveillance systems for Smart City applications. Existing grayscale mapping-based algorithms always suffer from over-enhancement of the background, noise amplification, and brightness distortion. To cope with these problems, an infrared image enhancement method based on adaptive histogram partition and brightness correction is proposed. First, the grayscale histogram is adaptively segmented into several sub-histograms by a locally weighted scatter plot smoothing algorithm and local minima examination. Then, the fore-and background sub-histograms are distinguished according to a proposed metric called grayscale density. The foreground sub-histograms are equalized using a local contrast weighted distribution for the purpose of enhancing the local details, while the background sub-histograms maintain the corresponding proportions of the whole dynamic range in order to avoid over-enhancement. Meanwhile, a visual correction factor considering the property of human vision is designed to reduce the effect of noise during the procedure of grayscale re-mapping. Lastly, particle swarm optimization is used to correct the mean brightness of the output by virtue of a reference image. Both qualitative and quantitative evaluations implemented on real infrared images demonstrate the superiority of our method when compared with other conventional methods.

1. Introduction

Infrared (IR) imaging has been extensively applied in Smart City applications [1,2], e.g., scene surveillance, goods sorting, and fire prevention due to its unique ability to receive IR rays (780 nm–300 $\mathsf{\mu }\mathrm{m}$) which cannot be perceived by human eyes. However, compared with visible (Vis) images, IR images do suffer from many intrinsic drawbacks, e.g., lack of color information, low contrast, blurred resolution, and visual disturbance caused by noise, which generate much inconvenience when attempting to recognize the target of interest from the background [3,4]. As a result, IR image enhancement is always a hot topic worthy of investigation in Smart City applications and plays a significant role in the pre-processing techniques of intelligent systems.

During the past decades, a number of image enhancement algorithms based on histogram equalization (HE) have been developed, the goal of which is to re-map the grayscales of the original image and obtain a new histogram with a uniform intensity distribution so as to improve the image contrast as much as possible. Among them, global histogram equalization (GHE) is the simplest one, in which a transfer function formulated by the cumulative density function (CDF) is utilized to re-calculate the grayscale histogram so that the overall histogram distribution is forced to be flattened out and accounts for a broader dynamic range, exaggerating the global contrast to the greatest extent [5]. To this end, grayscales with high probability distribution function (PDF) values (usually belonging to the background) will occupy the most intensity ranges and be dramatically enhanced, while those with low PDF values (usually belonging to target regions) will be greatly suppressed or even lost. That is to say, over-enhancement may easily occur in IR images with a large area of homogeneous regions [6]. In an even worse scenario, the noise existing in background may be hugely amplified due to over-enhancement.

Recently, some improved HE-based algorithms have also been investigated. Dynamic histogram equalization (DHE) proposed by Abdullah-Al-Wadud et al. [7] tends to preserve the details of the input. The raw grayscale histogram is partitioned based on local minima and specific grayscale ranges are assigned to each partition, after which HE is implemented on each sub-histogram [8]. Brightness preserving bi-histogram equalization (BBHE) [9] divides the original grayscale histogram into two parts using the mean brightness of the input and then the two sub-histograms are equalized independently. It has been proved that selecting the mean brightness as the separation threshold can preserve the general brightness before and after enhancement to a certain degree. Equal area dualistic sub-image histogram equalization (DSIHE) [10] increases the global contrast with the same strategy as BBHE, but the only difference is that the mean brightness is replaced by the median as the separation threshold. Relevant experiments further verify that DSIHE can achieve better performances than BBHE with respect to preserving the average information content of the input image. In addition, minimum mean brightness error bi-histogram equalization (MMBEBHE) [11] is another improvement of BBHE which yields the optimized mean brightness difference between the input and output. Recursive mean-separate histogram equalization (RMSHE) [12] and recursive sub-image histogram equalization (RSIHE) [13] are the recursive versions of BBHE and DSIHE. Although they provide a flexible way of monitoring the over-enhancement degree, the mean brightness is overly emphasized. Overall, the afore-introduced algorithms may have satisfactory enhancement results in Vis images, but side effects may happen when facing the variation of grayscale distribution in the histogram [14] and the strategy of brightness preservation is obviously unreasonable in IR image enhancement considering the fact that the input brightness itself is distorted.

With the aim of alleviating over-enhancement, Vicker et al. [15] proposed plateau histogram equalization (PHE), whose PDF is limited via a threshold. Furthermore, Song et al. [16] presented double plateau histogram equalization (DPHE) in which the upper threshold is designed to prevent the over-enhancement of background noise with typical grayscales, while the lower threshold is developed to protect details with fewer pixels from being merged [17]. Yet, the thresholds in PHE and DPHE need to be manually selected, which limits their usage in practice. As a result, adaptive plateau histogram equalization (APHE) [18], as well as adaptive double plateau histogram equalization (ADPHE) [17], were then developed. However, the robustness of ADPHE cannot be guaranteed because the thresholds are fixed for all the grayscales.

Most of the existing methods are highly dependent on mode estimation, which may easily lead to an over or under segmentation of the histogram. To deal with this significant problem, a non-parametric method using contrario theory-based automatic histogram mode estimation was proposed by Delon et al. [19]. In this method, a segmentation of the histogram without a priori assumptions in terms of the number or shape of its modes is presented, and the core idea is to test the simplest multi-modal law fitting the data. Specially speaking, an adequacy test named ‘meaningful rejections’ is developed, which is demonstrated to be both locally and globally superior. The algorithm is able to detect even very small modes when they are isolated, making it well adapted to document analysis. In addition, it is more robust to quantization noise owing to its statistical aspect.

To address the afore-mentioned problems in the existing methods, we focus on investigating an effective IR image enhancement algorithm based on adaptive histogram partition and brightness correction. First of all, the original grayscale histogram is smoothed by Gaussian filtering and a locally weighted scatter plot smoothing (LOWESS) algorithm [20] so that the spikes can be removed. After that, the histogram is partitioned into several intervals via local minima and whether a certain sub-histogram belongs to the foreground or background is judged by its grayscale density. For the foreground parts, a novel local contrast weighted distribution is presented to take the place of the conventional intensity distribution, which only takes the occurrence frequency of each gray level into account and each foreground sub-image is equalized using that distribution; for the background parts, all the corresponding intervals keep their proportions of the whole dynamic range. It is remarkable that, motivated by the characteristic of human eyes, a visual correction factor is proposed to adjust the range of each foreground histogram before equalization. Finally, a reference IR image with an appropriate brightness is employed to correct the mean brightness of the enhanced image by a particle swarm optimization (PSO) algorithm. Figure 1 shows a complete flow chart of our method.

2. Review of Global Histogram Equalization

Global histogram equalization (GHE) is the basic HE-based image enhancement method in computer vision and has been widely applied to practice. Its core is to build a transfer function based on PDF and re-map the raw grayscale distribution to a new uniform one.

Let us suppose that $\mathrm{I}=\{\mathrm{I}(\mathrm{x},\mathrm{y})\}$ represents the input IR image, where $\mathrm{I}(\mathrm{x},\mathrm{y})$ is the grayscale of each pixel whose spatial coordinate in $\mathrm{I}$ is $(\mathrm{x},\mathrm{y})$. Additionally, the IR image is digitized into L gray levels as $\{{\mathrm{I}}_0,{\mathrm{I}}_1,\dots ,{\mathrm{I}}_{\mathrm{L}-1}\}$ (for an 8-bit digital image, $\mathrm{L}=256$). So, it is apparent that $\forall \mathrm{I}(\mathrm{x},\mathrm{y})\in \{{\mathrm{I}}_0,{\mathrm{I}}_1,\dots ,{\mathrm{I}}_{\mathrm{L}-1}\}$. Next, the PDF of each grayscale can be defined as:

$$\mathrm{p}({\mathrm{I}}_{\mathrm{k}})=\frac{{\mathrm{n}}_{\mathrm{k}}}{\mathrm{n}}, \mathrm{k}=0,1,\dots ,\mathrm{L}-1,$$

(1)

where, ${\mathrm{I}}_{\mathrm{k}}$ refers to an arbitrary grayscale and $\mathrm{p}({\mathrm{I}}_{\mathrm{k}})$ is its corresponding PDF; ${\mathrm{n}}_{\mathrm{k}}$ denotes the number of pixels whose grayscales are ${\mathrm{I}}_{\mathrm{k}}$; and n is the total pixel number of I. Particularly, the graphic appearance of PDF is regarded as a histogram [14]. Further, the CDF of each grayscale is calculated as:

$$\mathrm{c}({\mathrm{I}}_{\mathrm{k}})={\displaystyle \sum _{\mathrm{t}=0}^{\mathrm{k}}\mathrm{p}({\mathrm{I}}_{\mathrm{t}}})={\displaystyle \sum _{\mathrm{t}=0}^{\mathrm{k}}\frac{{\mathrm{n}}_{\mathrm{t}}}{\mathrm{n}}},\mathrm{k}=0,1,\dots ,\mathrm{L}-1.$$

(2)

Obviously, ${\displaystyle \sum _{\mathrm{t}=0}^{\mathrm{L}-1}\mathrm{c}({\mathrm{I}}_{\mathrm{t}})=1}$. Based on Equation (2), the transfer function $\mathrm{T}({\mathrm{I}}_{\mathrm{k}})$ is constructed as:

$$\mathrm{T}({\mathrm{I}}_{\mathrm{k}})={\mathrm{I}}_0+({\mathrm{I}}_{\mathrm{L}-1}-{\mathrm{I}}_0)\cdot \mathrm{c}({\mathrm{I}}_{\mathrm{k}}).$$

(3)

Lastly, the enhanced image O is obtained as:

$$\mathrm{O}=\mathrm{T}(\mathrm{I})=\{\mathrm{T}(\mathrm{I}(\mathrm{x},\mathrm{y})|\forall \mathrm{I}(\mathrm{x},\mathrm{y})\in \{{\mathrm{I}}_0,{\mathrm{I}}_1,\dots ,{\mathrm{I}}_{\mathrm{L}-1}\})\}.$$

(4)

Note that:

$$\frac{\mathrm{d}\mathrm{c}({\mathrm{I}}_{\mathrm{k}})}{\mathrm{d}{\mathrm{I}}_{\mathrm{k}}}=\mathrm{c}({\mathrm{I}}_{\mathrm{k}})-\mathrm{c}({\mathrm{I}}_{\mathrm{k}-1})={\displaystyle \sum _{\mathrm{t}=0}^{\mathrm{k}}\mathrm{p}({\mathrm{I}}_{\mathrm{t}})-}{\displaystyle \sum _{\mathrm{t}=0}^{\mathrm{k}-1}\mathrm{p}({\mathrm{I}}_{\mathrm{t}})}=\mathrm{p}({\mathrm{I}}_{\mathrm{k}}).$$

(5)

From Equation (5), it is easy to find that grayscales with high PDFs can have a high gain, and account for the majority of the whole dynamic range. As a rule, those grayscales mainly belong to background pixels and are always over-merged by GHE, which is absolutely unexpected. On the contrary, the foreground, especially some small-sized targets, the grayscales of which have relatively low PDFs, are greatly suppressed after GHE, i.e., the details are thus lost.

3. Theory

In this section, the detailed theories related to the proposed algorithm are discussed at great length. First, Section 3.1 introduces the automatic way of partitioning the raw grayscale histogram. Next, the criterion of distinguishing between fore-and background sub-histograms is revealed in Section 3.2. After that, we clarify the different means of enhancing fore-and background sub-histograms, respectively, in Section 3.3. Finally, PSO algorithm optimizing the mean brightness with a reference image is depicted in Section 3.4.

3.1. Adaptive Segmentation of Input Grayscale Histogram

Considering the physical law in IR imaging that the temperature distribution of the target (foreground) or background is continuous and centers around a specific temperature, meaning that the grayscale distribution of the foreground should be unimodal (the condition of background is the same), we propose to judge whether a cluster of grayscales belongs to the foreground or background by means of partitioning the raw histogram into several intervals and recognizing each of them through a specific characteristic of grayscale distribution discussed in Section 3.2. In this way, we can use different strategies to enhance the fore-and background separately.

Inspired by previous works [21,22,23], the strategy of segmenting a histogram via local minima is adopted in this paper. Unlike the conventional methods, we focus on improving its robustness to those spikes existing in the histogram. Therefore, a data smoothing method called LOWESS is utilized as a pre-processing step.

We notice that there are large quantities of spikes, which can be seen as local outliers, in the input grayscale histograms, leading to too many peaks if the local minima are examined directly. Thus, a one-dimensional Gaussian filter is employed to roughly smooth the histogram at first:

$${\mathrm{p}}^{\prime }({\mathrm{I}}_{\mathrm{k}})={\displaystyle \sum _{\mathrm{j}=\mathrm{k}-{\mathrm{w}}_1}^{\mathrm{k}+{\mathrm{w}}_1}\mathrm{p}({\mathrm{I}}_{\mathrm{j}})}\cdot \mathsf{\kappa }(\mathrm{j}),\mathrm{k}={\mathrm{w}}_1,{\mathrm{w}}_1+1,\dots ,\mathrm{L}-{\mathrm{w}}_1-1,$$

(6)

where, ${\mathrm{p}}^{\prime }({\mathrm{I}}_{\mathrm{k}})$ refers to the smoothed PDF after Gaussian filtering; the size of the local window is $1*(2{\mathrm{w}}_1+1)$; and $\mathsf{\kappa }(\cdot )$ represents a one-dimensional Gaussian kernel function which is defined as:

$$\mathsf{\kappa }(\mathrm{x})=\frac{1}{\sqrt{2\mathsf{\pi }}{\mathsf{\sigma }}_1}{\mathrm{e}}^{-\frac{{(\mathrm{x}-{\mathrm{x}}_{\mathrm{o}})}^2}{2{\mathsf{\sigma }}_1{}^2}},$$

(7)

where, ${\mathsf{\sigma }}_1$ is a constant parameter deciding the weight of each neighboring PDF affecting the output and ${\mathrm{x}}_{\mathrm{o}}=\mathrm{k}$ denotes the index of the current grayscale ${\mathrm{I}}_{\mathrm{k}}$.

Then, a local smoother called LOWESS is used on the histogram after Gaussian filtering to further remove the spikes. The specific steps to nullify the minor peaks from ${\mathrm{p}}^{\prime }({\mathrm{I}}_{\mathrm{k}})$ are described as follows:

(1) Given a span $\mathsf{\zeta }$ with a fixed length $\mathrm{d}$, the regression weight of each neighboring PDF in $\mathsf{\zeta }$ is locally calculated using a 9-degree polynomial ${\mathsf{\phi }}_{\mathrm{j}}({\mathrm{p}}^{\prime }({\mathrm{I}}_{\mathrm{k}}))$ expressed as Equation (8).

$${\mathsf{\phi }}_{\mathrm{j}}({\mathrm{p}}^{\prime }({\mathrm{I}}_{\mathrm{k}}))={\left(1-{\left|\frac{{\mathrm{p}}^{\prime }({\mathrm{I}}_{\mathrm{k}})-{\mathrm{p}}^{\prime }({\mathrm{I}}_{\mathrm{j}})}{\mathrm{d}}\right|}^3\right)}^3,$$

(8)

where, ${\mathrm{p}}^{\prime }({\mathrm{I}}_{\mathrm{j}})$ denotes the neighboring PDF of ${\mathrm{p}}^{\prime }({\mathrm{I}}_{\mathrm{k}})$ in $\mathsf{\zeta }$ and ${\mathrm{p}}^{\prime }({\mathrm{I}}_{\mathrm{k}})$ is the current PDF to be smoothed.

(2) A polynomial unary regression estimation function is used to compute the smoothed PDF of each grayscale, and the coefficients of the locally linear regression are deduced by the weighted least square estimate of equations fitting an affine model to the probabilities. Please refer to the study made by Cleveland [20] for the detailed deduction. Here, the solution of the coefficients is directly given as:

$$\widehat{\mathsf{\theta }}=\left[\begin{array}{c}{\mathsf{\theta }}_0\\ {\mathsf{\theta }}_1\end{array}\right]={\left({\mathrm{X}}^{\mathrm{T}}\mathsf{\phi }\mathrm{X}\right)}^{-1}\left({\mathrm{X}}^{\mathrm{T}}\mathsf{\phi }\mathrm{Y}\right),$$

(9)

where, $\mathrm{X}\in {\mathrm{R}}^{\mathrm{d}\times 2}$ in which the first column values are equal to 1 and the second column values are written as a vector $({\mathrm{p}}^{\prime }({\mathrm{I}}_{\mathrm{k}})-{\mathrm{p}}^{\prime }({\mathrm{I}}_{\mathrm{j}}))$; $\mathsf{\phi }=\mathrm{diag}({\mathsf{\phi }}_1({\mathrm{p}}^{\prime }({\mathrm{I}}_{\mathrm{k}})),{\mathsf{\phi }}_2({\mathrm{p}}^{\prime }({\mathrm{I}}_{\mathrm{k}})),\dots ,{\mathsf{\phi }}_{\mathrm{d}}({\mathrm{p}}^{\prime }({\mathrm{I}}_{\mathrm{k}})))$ is a $\mathrm{d}\times \mathrm{d}$ sized diagonal matrix; and $\mathrm{Y}$ is the intended span of ${\mathrm{p}}^{\prime }({\mathrm{I}}_{\mathrm{k}})$. Figure 2 gives a simple schematic diagram of the above-discussed procedure of estimating.

According to Equation (9), a $2\times 1$ sized coefficient vector $\widehat{\mathsf{\theta }}$ of the locally linear regression model is yielded and the smoothed output ${\mathrm{p}}^{*}({\mathrm{I}}_{\mathrm{k}})$ of each PDF after Gaussian filtering ${\mathrm{p}}^{\prime }({\mathrm{I}}_{\mathrm{k}})$ is determined through Equation (10).

$${\mathrm{p}}^{*}({\mathrm{I}}_{\mathrm{k}})={\mathsf{\theta }}_1\cdot {\mathrm{p}}^{\prime }({\mathrm{I}}_{\mathrm{k}})+{\mathsf{\theta }}_0.$$

(10)

Figure 3 shows an example of smoothing the grayscale histograms produced by Gaussian filtering and LOWESS, respectively. It can be clearly seen that although the Gaussian-smoothed histogram in Figure 3b becomes smoother, there are still some fake peak and valley points (see the orange circle). Fortunately, in Figure 3c, LOWESS removes almost all the spikes effectively, while only the valid peaks and valleys remain. Thus, it proves that LOWESS is necessary in our work.

Based on the idea of searching for local minima, a $1\times {\mathrm{w}}_2$ sized sliding window is designed and all the local minima are extracted using the following criterion:

$$\mathrm{p}({\mathrm{I}}_{\mathrm{k}})=\{\begin{array}{c}\min \mathrm{imum},\mathrm{if}{\mathrm{p}}^{*}{(\mathrm{I}}_{\mathrm{k}})=\min ({\mathrm{p}}^{*}({\mathrm{I}}_{\mathrm{k}-({\mathrm{w}}_2-1)/2}),\dots ,{\mathrm{p}}^{*}({\mathrm{I}}_{\mathrm{k}}),\dots ,{\mathrm{p}}^{*}({\mathrm{I}}_{\mathrm{k}+({\mathrm{w}}_2+1)/2}))\hfill \\ \mathrm{non}-\min \mathrm{imum},\mathrm{others}\hfill \end{array},$$

(11)

where, $\min (\cdot )$ represents the minimal component in a set.

Besides, both the first and last non-zero components of the input grayscale histogram are also regarded as local minima. In this case, we assume that there is a total of $\mathsf{\varsigma }$ local minima in the histogram and they are denoted as $\{{\mathrm{m}}_1,{\mathrm{m}}_2,\dots ,{\mathrm{m}}_{\mathsf{\varsigma }}\}$. Hence, the raw histogram can be partitioned into $\mathsf{\varsigma }-1$ intervals, which are further written as $[{\mathrm{m}}_1,{\mathrm{m}}_2],[{\mathrm{m}}_2,{\mathrm{m}}_3],\dots ,[{\mathrm{m}}_{\mathsf{\varsigma }-1},{\mathrm{m}}_{\mathsf{\varsigma }}]$, and Figure 4 shows a more intuitive presentation of the segmentation result ($\mathsf{\varsigma }=7$).

3.2. Recognition of Fore-and Background Sub-Histograms

Before re-mapping the grayscales, we need to distinguish between fore- and background intervals and use different tools to deal with them. In this paper, we argue that whether a certain interval belongs to the foreground or background cannot depend on the PDF of a single grayscale, but the overall grayscale distribution of the interval should be considered. On account of the following two physical facts in IR imaging: (1) background pixels occupy a large percentage of the whole image whereas the object has an opposite condition, and both the intensity values of the foreground and background tend to cluster around a peak; (2) foreground intervals always present a flat shape whereas background intervals possess a narrow but sharp shape, a metric ${\Lambda }_{\mathrm{i}}$ measuring the grayscale density of each interval is proposed below:

$${\Lambda }_{\mathrm{i}}=\frac{{\mathrm{C}}_{\mathrm{i}}}{{\mathrm{N}}_{\mathrm{i}}},\mathrm{i}=1,2,\dots ,\mathrm{m}-1,$$

(12)

where, ${\mathrm{N}}_{\mathrm{i}}={\mathrm{m}}_{\mathrm{i}+1}-{\mathrm{m}}_{\mathrm{i}}+1$ refers to the quantity of grayscales in $[{\mathrm{m}}_{\mathrm{i}},{\mathrm{m}}_{\mathrm{i}+1}]$ and ${\mathrm{C}}_{\mathrm{i}}={\displaystyle \sum _{\mathrm{i}={\mathrm{m}}_{\mathrm{i}}}^{{\mathrm{m}}_{\mathrm{i}+1}}\mathrm{p}({\mathrm{I}}_{\mathrm{i}})}$ stands for the CDF of $[{\mathrm{m}}_{\mathrm{i}},{\mathrm{m}}_{\mathrm{i}+1}]$. Based on the definition of grayscale density, we further propose that a small ${\Lambda }_{\mathrm{i}}$ is related to the foreground and a large one corresponds to the background.

Next, an adaptive threshold ${\Lambda }^{*}$ to recognize the property of each interval is adaptively calculated via Ostu’s method [24], which is based on maximizing the inter-class variance of the dataset. Suppose that all the ${\Lambda }_{\mathrm{i}}$ compose a set $\Lambda =\{{\Lambda }_1,{\Lambda }_2,\dots ,{\Lambda }_{\mathsf{\varsigma }-1}\}$, and ${\Lambda }^{*}$ separates this set into two classes ${\Lambda }_{\mathsf{\alpha }}$ and ${\Lambda }_{\mathsf{\beta }}$. Ostu’s method exhaustively searches for the optimum ${\Lambda }^{*}$ as:

$${\Lambda }^{*}=\underset{{\Lambda }_{\min }<{\Lambda }^{*}<{\Lambda }_{\max }}{\arg \max }\{\mathsf{\sigma }({\Lambda }^{*})\},$$

(13)

where, ${\Lambda }_{\max }$ and ${\Lambda }_{\min }$ stand for the maximum and minimum of $\Lambda$, respectively; and,

$$\mathsf{\sigma }({\Lambda }^{*})={\mathrm{W}}_{\mathsf{\alpha }}\cdot {(\mathrm{E}({\Lambda }_{\mathsf{\alpha }})-{\mathsf{\mu }}_{\Lambda })}^2+{\mathrm{W}}_{\mathsf{\beta }}\cdot {(\mathrm{E}({\Lambda }_{\mathsf{\beta }})-{\mathsf{\mu }}_{\Lambda })}^2,$$

(14)

where, $\mathrm{E}({\Lambda }_{\mathsf{\alpha }})$ and $\mathrm{E}({\Lambda }_{\mathsf{\beta }})$ stand for the mean values of the two separated classes; ${\mathsf{\mu }}_{\Lambda }$ is the mean value of $\Lambda$; and ${\mathrm{W}}_{\mathsf{\alpha }}$ and ${\mathrm{W}}_{\mathsf{\beta }}$ represent the fractions indicating the component numbers of the two classes among the whole, respectively.

To this end, the grayscale histogram of input is segmented into fore- and background sub-histograms as the following criterion:

$$\{\begin{array}{c}{\Lambda }_{\mathrm{i}}\le {\Lambda }^{*},[{\mathrm{m}}_{\mathrm{i}},{\mathrm{m}}_{\mathrm{i}+1}]\in \mathrm{foreground}\\ {\Lambda }_{\mathrm{i}}>{\Lambda }^{*},[{\mathrm{m}}_{\mathrm{i}},{\mathrm{m}}_{\mathrm{i}+1}]\in \mathrm{background}\end{array},$$

(15)

Additionally, Figure 5 illustrates the above-introduced segmentation result.

3.3. Methods of Sub-Histogram Enhancement

In this section, some specific techniques are developed to enhance the partitioned sub-histograms respectively. Section 3.3.1 describes the details of allocating a mapped range for each interval; in Section 3.3.2, the fore-and background parts are enhanced with two different strategies.

3.3.1. Allocation of Re-Mapped Range

In order to overcome the adverse effect of over-enhancement, which commonly occurs in the background, we consider only equalizing the foreground sub-histograms while the modified intensity range of each background sub-histogram accounts for the same percentage in the raw histogram [25]. Under the circumstances, the gray levels of background will not be over-merged, but the target/background contrast is improved.

For an arbitrary background interval $[{\mathrm{m}}_{\mathrm{i}},{\mathrm{m}}_{\mathrm{i}+1}]$, its re-mapped range ${\mathrm{R}}_{\mathrm{i}}$ is computed as follows:

$${\mathrm{R}}_{\mathrm{i}}=\frac{{\mathrm{m}}_{\mathrm{i}+1}-{\mathrm{m}}_{\mathrm{i}}}{{\mathrm{I}}_{\max }-{\mathrm{I}}_{\min }}\cdot \mathrm{L},{\Lambda }_{\mathrm{i}}>{\Lambda }^{*},$$

(16)

where, ${\mathrm{I}}_{\max }$ and ${\mathrm{I}}_{\min }$ denote the maximal and minimal grayscales in the raw image, respectively.

The condition of foreground intervals is more complicated. First of all, the CDF is applied as a metric to decide the re-mapped range as:

$${\mathrm{R}}_{\mathrm{j}}=(\mathrm{L}-1-{\mathrm{m}}_1^{\prime }-{\displaystyle \sum _{\forall {\Lambda }_{\mathrm{s}}>{\Lambda }^{*}}{\mathrm{R}}_{\mathrm{s}}})\cdot \frac{{\mathrm{C}}_{\mathrm{j}}}{{\displaystyle \sum _{\forall {\Lambda }_{\mathrm{t}}\le {\Lambda }^{*}}{\mathrm{C}}_{\mathrm{t}}}},{\Lambda }_{\mathrm{j}}\le {\Lambda }^{*},$$

(17)

where, s and t represent the indexes of background and foreground intervals, respectively; ${\mathrm{C}}_{\mathrm{j}}={\displaystyle \sum _{\mathrm{k}={\mathrm{m}}_{\mathrm{j}}}^{{\mathrm{m}}_{\mathrm{j}+1}}\mathrm{p}({\mathrm{I}}_{\mathrm{k}})}$ is the CDF of a foreground interval $[{\mathrm{m}}_{\mathrm{j}},{\mathrm{m}}_{\mathrm{j}+1}]$; and ${\mathrm{m}}_1^{\prime }$ denotes the optimal minimal grayscale of the enhanced image, which will be discussed in Section 3.4.

In light of the property of human vision that noise in flat regions of an image will give rise to spurious or texture to the observer and the contrast sensitivity will decrease at sharp transitions [26], an exponential correction factor $\mathsf{\gamma }$ is designed to modify the re-mapped ranges of foreground intervals so as to reduce the effect of image noise as much as possible. Based on this consideration, we further notice that local entropy is an extensively-used physical metric in informatics that measures the richness of texture of an image block. Thus, we choose to use local entropy as a significant component to construct $\mathsf{\gamma }$. The definition of ${\mathsf{\gamma }}_{\mathrm{j}}\in [0,0.6]$ is given as:

$${\mathsf{\gamma }}_{\mathrm{j}}=\frac{{\Theta }_{\mathrm{j}}}{{\Theta }_{\max }}\times 0.6,$$

(18)

where, ${\Theta }_{\mathrm{j}}$ stands for the smoothness degree of $[{\mathrm{m}}_{\mathrm{j}},{\mathrm{m}}_{\mathrm{j}+1}]$ which is defined as Equation (19) and ${\Theta }_{\max }=\max \{{\Theta }_1,{\Theta }_2,\dots ,{\Theta }_{\mathsf{\varsigma }-1}\}$ is utilized for normalization.

$${\Theta }_{\mathrm{j}}=\frac{{\displaystyle \sum _{\mathrm{x}}{\displaystyle \sum _{\mathrm{y}}\mathrm{H}(\mathrm{x},\mathrm{y})}}}{{\mathrm{N}}_{\mathrm{j}}},(\mathrm{x},\mathrm{y})\in \{(\mathrm{x},\mathrm{y})|\mathrm{I}(\mathrm{x},\mathrm{y})\in [{\mathrm{m}}_{\mathrm{j}},{\mathrm{m}}_{\mathrm{j}+1}]\},$$

(19)

where, ${\mathrm{N}}_{\mathrm{j}}$ represents the number of pixels belonging to $[{\mathrm{m}}_{\mathrm{j}},{\mathrm{m}}_{\mathrm{j}+1}]$; and $\mathrm{H}(\mathrm{x},\mathrm{y})$ is the local entropy of pixel $(\mathrm{x},\mathrm{y})$ whose grayscale is within $[{\mathrm{m}}_{\mathrm{j}},{\mathrm{m}}_{\mathrm{j}+1}]$ and its formula is given as Equation (20).

$$\mathrm{H}(\mathrm{x},\mathrm{y})={\displaystyle \sum _{{\mathrm{W}}_{\mathrm{e}}(\mathrm{x},\mathrm{y})}{\displaystyle \sum _{\mathrm{t}=0}^{\mathrm{L}-1}-\mathrm{p}(\mathrm{t})\cdot {\log }_2}\mathrm{p}(\mathrm{t})\mathrm{i}\mathrm{f}\mathrm{p}(\mathrm{t})=0},{\log }_2\mathrm{p}(\mathrm{t})=0.$$

(20)

As is shown in Equation (18), the local entropy in regard to all the pixels in a ${\mathrm{w}}_3*{\mathrm{w}}_3$ sized local window ${\mathrm{W}}_{\mathrm{e}}(\mathrm{x},\mathrm{y})$ around (x,y) is counted. Obviously, a larger ${\Theta }_{\mathrm{j}}$ reveals a foreground interval with more complex textures on average.

As a result, Equation (17) can be re-written as:

$${\mathrm{R}}_{\mathrm{j}}=(\mathrm{L}-1-{\mathrm{m}}_1^{\prime }-{\displaystyle \sum _{\forall {\Lambda }_{\mathrm{s}}>{\Lambda }^{*}}{\mathrm{R}}_{\mathrm{s}}})\cdot \frac{{({\mathrm{C}}_{\mathrm{j}})}^{{\mathsf{\gamma }}_{\mathrm{j}}}}{{\displaystyle \sum _{\forall {\Lambda }_{\mathrm{t}}\le {\Lambda }^{*}}{({\mathrm{C}}_{\mathrm{t}})}^{{\mathsf{\gamma }}_{\mathrm{t}}}}},{\Lambda }_{\mathrm{j}}\le {\Lambda }^{*}.$$

(21)

3.3.2. Re-Mapping of Fore-and Background Histograms

Based on the principle that the background intervals maintain their proportions of the whole dynamic range without any other manipulations, the re-mapped result ${\mathrm{O}}^{\mathrm{b}}({\mathrm{I}}_{\mathrm{k}})$ of a background grayscale ${\mathrm{I}}_{\mathrm{k}}\in [{\mathrm{m}}_{\mathrm{i}},{\mathrm{m}}_{\mathrm{i}+1}]$ can be computed as Equation (22).

$${\mathrm{O}}^{\mathrm{b}}({\mathrm{I}}_{\mathrm{k}})={\displaystyle \sum _{\mathrm{j}=1}^{\mathrm{i}-1}{\mathrm{R}}_{\mathrm{j}}}+{\mathrm{R}}_{\mathrm{i}}\cdot \frac{{\mathrm{I}}_{\mathrm{k}}-{\mathrm{m}}_{\mathrm{i}}}{{\mathrm{m}}_{\mathrm{i}+1}-{\mathrm{m}}_{\mathrm{i}}}.$$

(22)

Then, let us further discuss the enhancement strategy for those foreground intervals, the local details of which need to be enhanced. In our study, a local contrast weighted distribution is developed and is written as:

$$\aleph ({\mathrm{I}}_{\mathrm{k}})=\frac{{\aleph }_{\mathrm{c}}({\mathrm{I}}_{\mathrm{k}})}{{\displaystyle \sum _{\mathrm{s}=0}^{\mathrm{L}-1}{\aleph }_{\mathrm{c}}({\mathrm{I}}_{\mathrm{s}})}},$$

(23)

where,

$${\aleph }_{\mathrm{c}}({\mathrm{I}}_{\mathrm{k}})={\displaystyle \sum _{\mathrm{x}=0}^{\mathrm{M}-1}{\displaystyle \sum _{\mathrm{y}=0}^{\mathrm{N}-1}\mathsf{\psi }(\mathsf{\omega }(\mathrm{x},\mathrm{y}))}}\cdot \mathsf{\delta }(\mathrm{I}(\mathrm{x},\mathrm{y}),{\mathrm{I}}_{\mathrm{k}}),$$

(24)

where, M and N are the width and height of the image, respectively; $\mathsf{\delta }(\mathrm{x},\mathrm{y})=\{\begin{array}{c}1,\mathrm{x}=\mathrm{y}\\ 0,\mathrm{others}\end{array}$ is the Kronecker delta function; and $\mathsf{\omega }(\mathrm{x},\mathrm{y})$ denotes the local contrast operator [27], which is defined as:

$$\mathsf{\omega }(\mathrm{x},\mathrm{y})=\frac{{\nabla }_{\mathrm{r}}(\mathrm{x},\mathrm{y})+{\nabla }_{\mathrm{l}}(\mathrm{x},\mathrm{y})+{\nabla }_{\mathrm{u}}(\mathrm{x},\mathrm{y})+{\nabla }_{\mathrm{d}}(\mathrm{x},\mathrm{y})}{4},$$

(25)

where,

$$\{\begin{array}{c}{\nabla }_{\mathrm{r}}(\mathrm{x},\mathrm{y})=\left|\mathrm{I}(\mathrm{x},\mathrm{y})-\mathrm{I}(\mathrm{x},\mathrm{y}+1)\right|\\ {\nabla }_{\mathrm{l}}(\mathrm{x},\mathrm{y})=\left|\mathrm{I}(\mathrm{x},\mathrm{y})-\mathrm{I}(\mathrm{x},\mathrm{y}-1)\right|\\ {\nabla }_{\mathrm{u}}(\mathrm{x},\mathrm{y})=\left|\mathrm{I}(\mathrm{x},\mathrm{y})-\mathrm{I}(\mathrm{x}+1,\mathrm{y})\right|\\ {\nabla }_{\mathrm{d}}(\mathrm{x},\mathrm{y})=\left|\mathrm{I}(\mathrm{x},\mathrm{y})-\mathrm{I}(\mathrm{x}-1,\mathrm{y})\right|\end{array},$$

(26)

As can be seen from Equations (25) and (26), $\mathsf{\omega }(\mathrm{x},\mathrm{y})$ indicates the four-direction gradient information of pixel $(\mathrm{x},\mathrm{y})$. Furthermore, $\mathsf{\psi }(\cdot )$ representing the weight function is formulated based on the sigmoid function as:

$$\mathsf{\psi }(\mathrm{x})=\mathsf{\alpha }\cdot \left(\mathrm{S}\mathrm{g}(\mathsf{\lambda }\cdot \mathrm{x}+\mathsf{\theta })+\mathsf{\xi }\right),$$

(27)

where, $\mathrm{S}\mathrm{g}(\mathrm{x})=\frac{1}{1+{\mathrm{e}}^{-\mathrm{x}}}$ is the standard sigmoid function whose curve is drawn as Figure 6a; and $\mathsf{\alpha }$, $\mathsf{\lambda }$, $\mathsf{\theta }$, and $\mathsf{\xi }$ denote the amplitude, scale, phase, and shift parameters, respectively.

Figure 6b reveals one of the desirable curve forms of $\mathsf{\psi }(\mathrm{x})$ that meets the following three requirements for our purpose: (1) a non-negative and monotonically increasing function; (2) $\underset{\mathrm{x}\to {\mathrm{x}}_{\min }}{\lim }\mathsf{\psi }(\mathrm{x})=0,{\mathrm{x}}_{\min }=\min (\mathsf{\omega })$; and (3) $\underset{\mathrm{x}\to {\mathrm{x}}_{\max }}{\lim }\mathsf{\psi }(\mathrm{x})=1,{\mathrm{x}}_{\max }=\max (\mathsf{\omega })$, and the exact values of the four parameters can be computed as:

$$\{\begin{array}{c}\mathsf{\alpha }=1\hfill \\ \mathsf{\lambda }=\frac{\ln \left(\frac{{\mathsf{\epsilon }}_2}{1-{\mathsf{\epsilon }}_2}\cdot \frac{1-{\mathsf{\epsilon }}_1}{{\mathsf{\epsilon }}_1}\right)}{\max (\mathsf{\omega })-\min (\mathsf{\omega })}\hfill \\ \mathsf{\theta }=\ln \left(\frac{{\mathsf{\epsilon }}_1}{1-{\mathsf{\epsilon }}_1}\right)-\frac{\ln \left(\frac{{\mathsf{\epsilon }}_2}{1-{\mathsf{\epsilon }}_2}\cdot \frac{1-{\mathsf{\epsilon }}_1}{{\mathsf{\epsilon }}_1}\right)}{\max (\mathsf{\omega })-\min (\mathsf{\omega })}\cdot \min (\mathsf{\omega })\hfill \\ \mathsf{\xi }=0\hfill \end{array},$$

(28)

where, ${\mathsf{\epsilon }}_1$ and ${\mathsf{\epsilon }}_2$ are two constant values that are approach 0 and 1, respectively. An example of calculating the local contrast weighted distribution via the above-introduced method is given in Figure 7.

On the basis of the presented local contrast weighted distribution $\aleph ({\mathrm{I}}_{\mathrm{k}})\in [0,1]$, the grayscale distribution of the foreground interval is substituted and a foreground grayscale ${\mathrm{I}}_{\mathrm{k}}\in [{\mathrm{m}}_{\mathrm{j}},{\mathrm{m}}_{\mathrm{j}+1}]$ is re-mapped using the approach expressed as:

$${\mathrm{O}}^{\mathrm{f}}({\mathrm{I}}_{\mathrm{k}})={\displaystyle \sum _{\mathrm{t}=1}^{\mathrm{j}-1}{\mathrm{R}}_{\mathrm{t}}}+{\mathrm{R}}_{\mathrm{j}}\cdot \frac{\aleph ({\mathrm{I}}_{\mathrm{k}})}{{\mathrm{C}}_{\mathrm{j}}^{\aleph }},$$

(29)

where, ${\mathrm{C}}_{\mathrm{j}}^{\aleph }={\displaystyle \sum _{\mathrm{k}={\mathrm{m}}_{\mathrm{j}}}^{{\mathrm{m}}_{\mathrm{j}+1}}\aleph ({\mathrm{I}}_{\mathrm{k}})}$ denotes the cumulative local contrast weighted distribution of $[{\mathrm{m}}_{\mathrm{j}},{\mathrm{m}}_{\mathrm{j}+1}]$. Here, a schematic diagram of the grayscale re-mapping is also drawn in Figure 8.

Thus, the final enhanced image O can be obtained as:

$$\mathrm{O}=\{\mathrm{O}(\mathrm{i},\mathrm{j})\}={\mathrm{O}}^{\mathrm{f}}\cup {\mathrm{O}}^{\mathrm{b}}.$$

(30)

However, one point that should be noticed is that the minimal grayscale after re-mapping ${\mathrm{m}}_1^{\prime }$ is uncertain and it decides the general brightness of the output image directly. The optimization of this parameter is discussed in Section 3.4.

3.4. Brightness Correction via PSO Algorithm

As is mentioned above, the minimal grayscale ${\mathrm{m}}_1^{\prime }$ affects the overall brightness of the enhanced IR image. In normal conditions, ${\mathrm{m}}_1^{\prime }$ is set as 0 by default, which results in the phenomenon that the enhanced image is over-dark. However, we cannot simply use the strategies of brightness preservation proposed in the previous works [9,10,11,12,13], which aim to make the mean brightness close to the input as much as possible, because the brightness of the input itself is distorted for IR images. On account of this problem, we consider using the PSO algorithm to optimize the exact value of ${\mathrm{m}}_1^{\prime }$ by virtue of a reference image with a suitable brightness.

PSO is a population-based stochastic optimization algorithm [28] and it utilizes a certain number of particles forming a swarm motion in the search space to estimate the optimum solution of the objective function F. All the particles fly through the search space depending on the best solution of each particle achieved so far and the best solution tracked by any particle among all generations of the swarm, which are named $\mathrm{pbest}$ and $\mathrm{gbest}$, respectively [29].

3.4.1. Particle Initialization

In this step, N particles are first generated. The initial position of each particle can be regarded as random sampling, and the position of each particle is a potential value of ${\mathrm{m}}_1^{\prime }$. Assuming that the sum of the re-mapping ranges of background intervals calculated in Section 3.3.1 is denoted as $\Re ={\displaystyle \sum _{\forall {\Lambda }_{\mathrm{s}}>{\Lambda }^{*}}{\mathrm{R}}_{\mathrm{s}}}$, we propose that the initial position of each particle should be greater than 0 and smaller than $\mathrm{L}-1-\Re$. Based on this point, uniform distribution-based random sampling is implemented in the range $[0,\mathrm{L}-1-\Re ]$ to assign the initial value for each particle.

Here, an example of particle initialization for Figure 9 is shown in

Table 1. Examples of particle initialization for Figure 9.

‘building’	Particle index	1	2	3	4	5	6	7	8	9	10
	Initial position	1	42	102	125	41	114	50	162	7	105
‘room’	Particle index	1	2	3	4	5	6	7	8	9	10
	Initial position	94	17	114	152	106	125	80	104	86	151

. $\Re =90$ for Figure 9a and $\Re =99$ for Figure 9b, so the sampling ranges are $[0,165]$ and $[0,156]$, respectively. Note that all the initial positions are sampled as integers in our study.

3.4.2. Objective Function

In our study, ${\mathrm{m}}_1^{\prime }$ to be optimized composes a one-dimensional search space for PSO and the objective function is defined as the absolute mean brightness error [30] of the enhanced image and reference image. Equation (29) gives the formulation of the objective function $\mathrm{F}$ as:

$$\mathrm{F}=\left|{\mathrm{M}}_{\mathrm{o}}({\mathrm{m}}_1^{\prime })-{\mathrm{M}}_{\mathrm{R}}\right|$$

(31)

where, ${\mathrm{M}}_{\mathrm{O}}({\mathrm{m}}_1^{\prime })$ and ${\mathrm{M}}_{\mathrm{R}}$ stand for the mean of output image O and reference image R, respectively. Unlike the offline learning strategy, the reference image R is not used for the learning purpose. Here, only a reference image, the mean brightness of which is comfortable for the human vision system, is needed, i.e., its mean brightness is neither too dark nor too bright. It is utilized to generate a standard mean intensity value ${\mathrm{M}}_{\mathrm{R}}$ which is regarded as a standard value to optimize the objective function value F in the following PSO optimization. By this means, the optimal value of ${\mathrm{m}}_1^{\prime }$ can be computed. Note that ${\mathrm{M}}_{\mathrm{O}}({\mathrm{m}}_1^{\prime })$ can be calculated using the method discussed from Section 3.1, Section 3.2 and Section 3.3, once ${\mathrm{m}}_1^{\prime }$ is generated by PSO.

3.4.3. Implementation of PSO Iteration

Every particle tries to change its position via information including the current position, current velocity, the respective distance between the current position and $\mathrm{pbest}$, and the distance between the current position and the $\mathrm{gbest}$ [31]. The velocity of each particle is suitably updated as Equation (30).

$${\mathrm{V}}_{\mathrm{i}}^{\mathrm{t}+1}=\mathsf{\omega }\cdot {\mathrm{V}}_{\mathrm{i}}^{\mathrm{t}}+{\mathrm{c}}_1\cdot \mathrm{rand}()\cdot ({\mathrm{pbest}}_{\mathrm{i}}-{\mathrm{P}}_{\mathrm{i}}^{\mathrm{t}})+{\mathrm{c}}_2\cdot \mathrm{rand}()\cdot ({\mathrm{gbest}}_{\mathrm{i}}-{\mathrm{P}}_{\mathrm{i}}^{\mathrm{t}}),$$

(32)

where, ${\mathrm{V}}_{\mathrm{i}}^{\mathrm{t}}$ denotes the velocity of the i-th particle at the t-th iteration; ${\mathrm{c}}_1$ and ${\mathrm{c}}_2$ are learning factors; $\mathrm{rand}()\in [0,1]$ represents a uniformly distributed random number; ${\mathrm{P}}_{\mathrm{i}}^{\mathrm{t}}$ is the position of the i-th particle at th t-th iteration; and ${\mathrm{pbest}}_{\mathrm{i}}$ represents the current best of the i-th particle (the minimal F in the current iteration) while ${\mathrm{gbest}}_{\mathrm{i}}$ represents the global best of the swarm (the minimal F in all the iterations that have been implemented so far). Besides, $\mathsf{\omega }$ is the inertia weight which is set to be in the range $[{\mathsf{\omega }}_{\min },{\mathsf{\omega }}_{\max }]$ and Equation (33) gives its formulation in detail.

$$\mathsf{\omega }=\frac{{\mathsf{\omega }}_{\max }-({\mathsf{\omega }}_{\max }-{\mathsf{\omega }}_{\min })\cdot \mathrm{t}}{{\mathrm{t}}_{\max }},$$

(33)

where, t is the current iteration and ${\mathrm{t}}_{\max }$ is the maximal iteration.

Based on the afore-discussed modified velocity ${\mathrm{V}}_{\mathrm{i}}^{\mathrm{t}+1}$, the position of a particle ${\mathrm{P}}_{\mathrm{i}}^{\mathrm{t}+1}$ can be thus updated as:

$${\mathrm{P}}_{\mathrm{i}}^{\mathrm{t}+1}={\mathrm{P}}_{\mathrm{i}}^{\mathrm{t}}+{\mathrm{V}}_{\mathrm{i}}^{\mathrm{t}+1}.$$

(34)

Then, the optimal ${\mathrm{m}}_1^{*}$ is achieved by iteration until convergence, and the complete procedure is concluded in Algorithm 1. Finally, the enhanced image is outputted using the procedures introduced in Section 3.1, Section 3.2 and Section 3.3.

Algorithm 1: Optimization of${\mathrm{m}}_1^{\prime }$using PSO

1. Input: an IR image $\mathrm{I}$ to be enhanced and the mean of a reference image ${\mathrm{M}}_{\mathrm{R}}$.

2. Create N particles $\{{\mathrm{P}}_1,{\mathrm{P}}_2,\dots ,{\mathrm{P}}_{\mathrm{N}}\}$, where ${\mathrm{P}}_{\mathrm{i}}\in {\mathrm{R}}^{1\times 1}$ representing an ${\mathrm{m}}_1^{\prime }$;

3. for each particle i = 1 to N do

4. Initialize the position of each particle ${\mathrm{P}}_{\mathrm{i}}^1$ and its corresponding velocity ${\mathrm{V}}_{\mathrm{i}}^1$;

5. end for

6. for t = 1 : ${\mathrm{t}}_{\max }$ do

7. for each particle

8.  Produce the enhanced image ${\mathrm{O}}_{\mathrm{i}}^{\mathrm{t}}$ using methods introduced in Section 3.1, Section 3.2 and Section 3.3 based on its position value ${\mathrm{P}}_{\mathrm{i}}^{\mathrm{t}}$;

9.  Calculate the objective function ${\mathrm{F}}_{\mathrm{i}}^{\mathrm{t}}$ using Equation (29);

10.  if ${\mathrm{F}}_{\mathrm{i}}^{\mathrm{t}}<\mathrm{F}({\mathrm{pbest}}_{\mathrm{i}})$

11.   ${\mathrm{pbest}}_{\mathrm{i}}={\mathrm{P}}_{\mathrm{i}}^{\mathrm{t}}$;

12.  end if

13.   if ${\mathrm{F}}_{\mathrm{i}}^{\mathrm{t}}<\mathrm{F}({\mathrm{gbest}}_{\mathrm{i}})$

14.   ${\mathrm{gbest}}_{\mathrm{i}}={\mathrm{P}}_{\mathrm{i}}^{\mathrm{t}}$;

15.   end if

16. end for

17. for each particle i = 1 : N do

18.  Update the velocity ${\mathrm{V}}_{\mathrm{i}}^{\mathrm{t}+1}$ using Equation (32);

19.  Update the position ${\mathrm{P}}_{\mathrm{i}}^{\mathrm{t}+1}$ using Equation (34);

20. end for

21. end for

22. Output: the optimal ${\mathrm{m}}_1^{*}$.

4. Experiment and Discussion

In this section, a number of real IR images are applied to make comparative experiments and they can be downloaded from a publicly available database [32]. The reference image with a suitable brightness used in our method is shown in Figure 10 and its mean intensity is 99. All the test IR images whose wavelengths on the electromagnetic spectrum are from 4 to 12 $\mathsf{\mu }\mathrm{m}$ [33] are 8-bit and their sizes are $384\times 288$. Figure 11 presents the eight test images utilized, which are named as ‘residence’, ‘car’, ‘guardrail’, ‘building’, ‘room’, ‘cabinet’, ‘leg’, and ‘computer’, respectively. The former fours are outdoor scenes, while the latter fours are indoor scenes. In our experiments, seven conventional algorithms: GHE, BBHE, MMBEBHE, DSIHE, RMSHE, RSIHE, and ADPHE, introduced in Section 1, are used to make both qualitative and quantitative comparisons. Among them, the recursion level of RMSHE is set as the default value $\mathrm{r}=2$.

Below, we firstly discuss the parameter setting in Section 4.1; the processing results as well as some related discussions are revealed in Section 4.2; in Section 4.3, several quantitative metrics are employed to further demonstrate the superiority of our method. Furthermore, all the simulations are implemented by Matlab 2012a software on a PC with a 2.60 GHz INTEL CPU and 4 GB memory. The codes of this work can be downloaded from our github website [34].

4.1. Parameter Setting

In this section, all the parameters utilized, as well as the corresponding settings, are listed in

Table 2. Parameter settings.

Parameter	Meaning	Default Value
${\mathrm{w}}_1$	Radius of Gaussian filtering window	4
${\mathsf{\sigma }}_1$	Variance of Gaussian filter	0.7
d	Length of LOWESS span	11
${\mathrm{w}}_2$	Length of sliding window for local minima examination	9
${\mathrm{w}}_3$	Length of local entropy window	7
${\mathsf{\epsilon }}_1$	Lower bound of the weight function $\aleph (\mathrm{x})$	0.0001
${\mathsf{\epsilon }}_2$	Upper bound of the weight function $\aleph (\mathrm{x})$	0.9999
${\mathrm{t}}_{\max }$	Maximal iteration of PSO	10
$\mathrm{N}$	Particle number of PSO	10
${\mathrm{c}}_1$	Learning factor of PSO	0.5
${\mathrm{c}}_2$	Learning factor of PSO	0.5
${\mathsf{\omega }}_{\max }$	Maximal inertia weight of PSO	0.9
${\mathsf{\omega }}_{\min }$	Minimal inertia weight of PSO	0.1

. The default values of ${\mathrm{c}}_1$, ${\mathrm{c}}_2$, ${\mathsf{\omega }}_{\max }$, and ${\mathsf{\omega }}_{\min }$ are set according to [31], and the others are set empirically. Note that all these values are fixed in all the simulation experiments.

Specifically speaking, ${\mathrm{w}}_1$ decides the length of the Gaussian filtering window and ${\mathsf{\sigma }}_1$ affects the weight of each neighboring PDF. The suitable values of these two parameters need to ensure that the remarkable spikes in the raw histogram can be filtered while the local feature of normal PDFs can be preserved. In the LOWESS algorithm, d adjacent PDFs are utilized to codetermine the final smoothed output of a centered PDF processed by Gaussian filtering. For the purpose of maintaining the local shape of the raw histogram, d cannot be too large. With regards to ${\mathrm{w}}_2$, it should be neither too large nor too small so that those spikes which cannot be totally smoothed by Gaussian filtering and LOWESS are prevented from being recognized as local minima and only those outstanding valley points can be picked out. We argue that ${\mathrm{w}}_3$ should be relatively small, because a large ${\mathrm{w}}_3$ would not only increase the computational cost, but also make the local entropy calculated not well represent the property of the centered pixel. Due to the fact that the weight function $\aleph (\mathrm{x})$ cannot exactly reach 0 and 1, ${\mathsf{\epsilon }}_1$ and ${\mathsf{\epsilon }}_2$ are used to represent the practical lower and upper bounds, respectively. Thus, their values should approach to 0 and 1 as much as possible. In order to limit the computational load and get a satisfactory convergence simultaneously, the maximal iteration ${\mathrm{t}}_{\max }$ and the particle number $\mathrm{N}$ are both set as 10. Lastly, based on previous work on the PSO algorithm [31], the two learning factors ${\mathrm{c}}_1$ and ${\mathrm{c}}_2$ should be in the range $[0,4]$, and ${\mathsf{\omega }}_{\max }$, ${\mathsf{\omega }}_{\min }$ are commonly set as 0.9 and 0.1, respectively.

4.2. Qualitative Comparison

First of all, the enhanced images of all the eight algorithms are given in Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19. As is shown in Figure 11, all the IR images to be enhanced suffer from low global contrast and blurred details, and image noise is particularly serious in Figure 11e–h, which are indoor scenes.

For the ‘residence’ scene, it can be observed from Figure 12 that GHE, BBHE, and MMBEBHE have strong abilities of grayscale equalization, leading to the exaggerated fore- and background contrast. The wall and the roof are thus over-enhanced and the noise on the ground is also obviously amplified. On the other hand, DISHE, RMSHE, and RSIHE maintain the input brightness to a certain extent, whereas some grayscales are still over merged. The details, e.g., the windows and doors, are enhanced by ADPHE, but the visual contrast is not clearly improved, which is far from our expectation. Fortunately, our proposed algorithm achieves a relatively better visual performance which avoids the over-enhancement on the wall successfully.

With respect to the ‘car’ image shown in Figure 11b, the thermal radiations of the wheel and snow cause severe disturbances for most of the traditional methods. The grayscales belonging to the wheel are overly improved and over-merged in GHE, BBHE, MMBEBHE, DISHE, and ADPHE, while those belonging to the telegraph pole and wall are mistakenly re-mapped by RMSHE and RSIHE, but all these adverse phenomena do not happen in our method.

The grass and guardrail pixels of the ‘guardrail’ scene are exceedingly suppressed by GHE, BBHE, and MMBEBHE, which results in an unnatural visual performance. Besides, artifacts are seriously generated in the wall region at the top left corner, especially in DISHE, RMSHE, RSIHE, and ADPHE (see Figure 14d–g). For our algorithm, we argue that although the global contrast of Figure 14h is lower than Figure 14a–c,g, the over-enhancement is completely overcome and the overall clarity is satisfactory.

The condition of the ‘building’ scene is quite similar to ‘residence’, and the transition area around the skyline becomes visually abrupt due to the fact that the conventional algorithms treat fore- and background in the same way. In addition, the image noise is comparatively obvious in Figure 15c produced by MMBEBHE and the textures of the block wall are lost dramatically in Figure 15g produced by ADPHE. It is worth noting that the enhanced result using RSIHE is much more blurred when compared with DISHE and RMSHE, even though they can all preserve the input brightness quite well, indicating that more detailed information is missing in RSIHE. When considering the performance of our method, we find that not only the troublesome noise is suppressed, but also the fore- and background contrast can meet the observational requirement well.

The ‘room’ image displayed in Figure 11e is changeling on account of the striking noise. On the one hand, the grayscales belonging to the box at the bottom right corner are overly enhanced, except in our method. What is more, the noise distributing around the door is extremely magnified by GHE, BBHE, and MMBEBHE. We also notice that the small heat source near the chairs is over-merged by ADPHE, despite the fact that the noise is not amplified.

The cabinet in Figure 11f is fuzzy and all the comparing algorithms make attempts to improve its contrast as much as possible. We consider this kind of strategy as incorrect since it decreases the clarity in spite of the improvement of contrast. In other words, the balance between clarity and contrast is the key point in IR image enhancement. As is shown in Figure 17h, our presented method not only preserves the texture of the floor, shelf, and power lines, but also increases the global contrast to a certain degree.

The ‘leg’ image has a relatively higher target/background contrast when it is compared with other test images. One of the most significant tasks in enhancing this image is to enrich the texture existing on the leg and to keep the local contrast of background. Among all the processed results, severe over-enhancement of grayscales belonging to the chair leg and the wall occurs in GHE, MMBEBHE, and ADPHE, even though ADPHE adaptively selects the upper threshold to suppress noise to some extent. It should be noted that although BBHE and DISHE achieve satisfactory global contrasts, one of the legs of the chair located in the middle of the image disappears and the noise is magnified more seriously than RMSHE, RSIHE, ADPHE, and our method.

Finally, the ‘computer’ image is the most challenging one among all the eight images due to the remarkable noise and quite weak contrast. It is clear that only our method succeeds in improving the image contrast while suppressing the noise. Owing to the usage of local contrast weighted distribution, the heat source on the table becomes more outstanding, while it totally vanishes in ADPHE. On the other hand, our strategy of maintaining the raw proportion of background interval in the input grayscale histogram makes it a success in simultaneously eliminating background noise and improving the global contrast.

4.3. Quantitative Comparison

4.3.1. Comparison of Processed Result

To further verify the superiority of the presented algorithm, several extensively used metrics are applied to make quantified comparisons in this section.

First of all, the definitions of all the metrics utilized in this paper are introduced:

(1) Linear index of fuzziness $\mathsf{\eta }$ [35] is widely applied to analyze the performance of image enhancement, which is defined as:

$$\mathsf{\eta }(\mathrm{O})=\frac{2}{\mathrm{M}\times \mathrm{N}}{\displaystyle \sum _{\mathrm{x}=0}^{\mathrm{M}-1}{\displaystyle \sum _{\mathrm{y}=0}^{\mathrm{N}-1}\min (\mathsf{\Psi }(\mathrm{x},\mathrm{y}),1-\mathsf{\Psi }(\mathrm{x},\mathrm{y}))}},$$

(35)

where, O represents the enhanced image and its width and height are M and N, respectively. and,

$$\mathsf{\Psi }(\mathrm{x},\mathrm{y})=\sin \left[\frac{\mathsf{\pi }}{2}\left(1-\frac{\mathrm{O}(\mathrm{x},\mathrm{y})}{\max (\mathrm{O})}\right)\right],$$

(36)

where, max(O) represents the maximal grayscale of O. A smaller value of $\mathsf{\rho }$ indicates a better performance of the enhanced image with fewer clutters [36].

(2) Peak signal-to-noise ratio (PSNR) [37,38,39] is commonly used for evaluating image quality and measures how much the enhanced image has degraded when referred to the input image in the area of image enhancement [38]. Here, the formula of PSNR is given in Equation (37).

$$\mathrm{PSNR}(\mathrm{O})=20\times \lg \left(\frac{\mathrm{MAX}}{\sqrt{\mathrm{MSE}(\mathrm{O},\mathrm{I})}}\right),$$

(37)

where, $\mathrm{MAX}$ denotes the maximal gray level and $\mathrm{MAX}=255$ for 8-bit digital images; MSE(O,I) is the mean square error between the output O and input I, which is defined as:

$$\mathrm{MSE}(\mathrm{O},\mathrm{I})=\frac{1}{\mathrm{M}\times \mathrm{N}}{\displaystyle \sum _{\mathrm{x}=0}^{\mathrm{M}-1}{\displaystyle \sum _{\mathrm{y}=0}^{\mathrm{N}-1}{(\mathrm{O}(\mathrm{x},\mathrm{y})-\mathrm{I}(\mathrm{x},\mathrm{y}))}^2}}.$$

(38)

Besides, a larger PSNR refers to a higher image quality of the output.

(3) Image definition $\mathsf{\rho }$ [40] is a comprehensive index integrating $\mathsf{\eta }$ and PSNR and is always used to reflect the overall definition of an enhanced image. Below, its specific definition is given as:

$$\mathsf{\rho }(\mathrm{O})=\frac{\mathrm{PSNR}(\mathrm{O})}{\mathsf{\eta }(\mathrm{O})}.$$

(39)

Ideally, a qualified result should not only improve the global contrast, but also has less degradation. Thus, for a high-quality enhanced image, a small $\mathsf{\rho }$ is expected.

(4) Roughness $\mathsf{\Omega }$ [41] is used as a metric for evaluating the performance of an algorithm for noise reduction in many IR imaging applications. Here, we employ it in the following experiments for measuring the effect of noise suppression of the enhanced image. Further, the calculation of $\mathsf{\Omega }$ is introduced below:

$$\mathsf{\Omega }(\mathrm{O})=\frac{{\Vert {\mathrm{h}}_1\ast \mathrm{O}\Vert }_1+{\Vert {\mathrm{h}}_2\ast \mathrm{O}\Vert }_1}{{\Vert \mathrm{O}\Vert }_1},$$

(40)

where, ${\mathrm{h}}_1=[1,-1]$ and ${\mathrm{h}}_2={[1,-1]}^{\mathrm{T}}$ stand for the horizonal and vertical difference filters, respectively; $\ast$ is a convolution operator; and ${\Vert \cdot \Vert }_1$ stands for the ${\mathrm{L}}_1$ norm. Obviously, a smaller $\mathsf{\Omega }$ is related to a smoother result with less residual noise.

(5) Discrete entropy (DE) [42] characterizing the information amount contained in an enhanced image is employed to measure the degree of over-enhancement in our experiment. Additionally, DE is a globally statistical index which is defined as:

$$\mathrm{D}\mathrm{E}(\mathrm{O})={\displaystyle \sum _{\mathrm{s}=0}^{\mathrm{L}-1}-\mathrm{p}(\mathrm{s})\cdot {\log }_2}\mathrm{p}(\mathrm{s})\mathrm{if}\mathrm{p}(\mathrm{s})=0,{\log }_2\mathrm{p}(\mathrm{s})=0.$$

(41)

A larger DE means fewer gray levels are merged, leading to a clearer visual performance. Note that the essence of DE is different from the one of local entropy H defined in Equation (20), although their formulas are seemingly the same. DE is a globally counted index, but H is calculated in a local window.

(6) Logarithmic Michelson Contrast Measure (AME) [43] denotes a measure of local contrast, and Equation (42) gives its definition as:

$$\mathrm{AME}(\mathrm{O})=-\frac{1}{{\mathrm{l}}_1\times {\mathrm{l}}_2}{\displaystyle \sum _{\mathrm{i}=1}^{{\mathrm{l}}_1}{\displaystyle \sum _{\mathrm{j}=1}^{{\mathrm{l}}_2}20\cdot \ln \left(\frac{{\mathrm{O}}_{\max }^{\mathrm{i},\mathrm{j}}-{\mathrm{O}}_{\min }^{\mathrm{i},\mathrm{j}}+\mathsf{\epsilon }}{{\mathrm{O}}_{\max }^{\mathrm{i},\mathrm{j}}+{\mathrm{O}}_{\min }^{\mathrm{i},\mathrm{j}}+\mathsf{\epsilon }}\right)}},$$

(42)

where, the output O is segmented into ${\mathrm{l}}_1\times {\mathrm{l}}_2$ blocks; ${\mathrm{O}}_{\max }^{\mathrm{i},\mathrm{j}}$ and ${\mathrm{O}}_{\min }^{\mathrm{i},\mathrm{j}}$ denote the maximal and minimal grayness of the block, respectively; and $\mathsf{\epsilon }={10}^{-5}$ is a constant for preventing the invalid values. AME aims to use the relationship between the spread and the sum of the two intensity values in each block, and a smaller AME means a better performance [33].

Table 3. Quantitative comparisons of $\mathsf{\eta }$ for test images, and bold values indicate the top two results.

	GHE	BBHE	MMBEBHE	DISHE	RMSHE	RSIHE	ADPHE	Ours
‘residence’	0.1974	0.1993	0.2094	0.2508	0.2546	0.2687	0.3188	0.1427
‘car’	0.1973	0.1871	0.2069	0.2446	0.2441	0.2589	0.2991	0.1389
‘guardrail’	0.1980	0.1986	0.2052	0.2500	0.2492	0.2642	0.2863	0.1196
‘building’	0.1979	0.1935	0.2096	0.2484	0.2494	0.2670	0.3113	0.1479
‘chair’	0.1978	0.1845	0.2033	0.2496	0.2405	0.2592	0.2718	0.1904
‘table’	0.1965	0.2132	0.2041	0.2530	0.2557	0.2684	0.2768	0.1484
‘leg’	0.1973	0.1537	0.2116	0.2336	0.2279	0.2496	0.2235	0.1204
‘computer’	0.2009	0.2082	0.2105	0.2412	0.2514	0.2766	0.2767	0.2228
Ave.	0.1978	0.1922	0.2075	0.2464	0.2466	0.2640	0.2625	0.1542

,

Table 4. Quantitative comparisons of PSNR for test images, and bold values indicate the top two results.

	GHE	BBHE	MMBEBHE	DISHE	RMSHE	RSIHE	ADPHE	Ours
‘residence’	11.8826	12.0603	12.1058	15.5140	14.2222	17.2296	18.8433	18.5197
‘car’	13.1067	13.4511	13.3919	16.8303	16.8147	20.5077	15.3592	16.9767
‘guardrail’	13.0533	13.2453	13.3280	16.6305	16.1224	19.4358	13.1659	16.7550
‘building’	12.2800	12.6611	12.7154	16.0500	15.3751	18.5320	16.8048	16.8685
‘chair’	13.1027	13.2009	13.2610	16.8133	16.1347	18.2692	12.5104	20.5891
‘table’	12.4319	12.5129	12.6715	16.1413	14.4491	17.9370	12.0847	16.5403
‘leg’	13.3853	14.2202	13.6672	19.2938	24.1548	24.4557	25.2335	24.6008
‘computer’	11.4300	11.4154	11.4812	14.9679	13.0504	16.4078	27.6319	21.4629
Ave.	12.5840	12.8459	12.82775	16.5301	16.2904	19.0968	17.7042	19.0391

,

Table 5. Quantitative comparisons of $\mathsf{\rho }$ for test images, and bold values indicate the top two results.

	GHE	BBHE	MMBEBHE	DISHE	RMSHE	RSIHE	ADPHE	Ours
‘residence’	60.1955	60.5132	57.8118	61.8580	55.8609	64.1220	59.1070	129.7807
‘car’	66.4303	71.8925	64.7264	68.8074	68.88447	79.2108	51.3513	122.2224
‘guardrail’	65.9257	66.6933	64.9512	66.5220	64.6966	73.5647	45.9863	140.0919
‘building’	62.0515	65.4320	60.6650	64.6135	61.6483	69.4082	53.9826	114.0534
‘chair’	66.2421	71.5495	65.2287	67.3609	67.0881	70.4830	46.0279	108.1360
‘table’	63.2666	58.6909	62.0847	63.7996	56.5080	66.8293	43.6585	111.4575
‘leg’	67.8423	92.5191	64.5897	82.5933	105.9885	97.9795	112.9015	204.3256
‘computer’	56.8939	54.8290	54.5425	62.0559	51.9108	59.3195	99.8623	96.3326
Ave.	63.6059	67.7649	61.8250	67.2013	66.5732	72.6146	64.1096	128.3000

,

Table 6. Quantitative comparisons of $\mathsf{\Omega }$ for test images, and bold values indicate the top two results.

	GHE	BBHE	MMBEBHE	DISHE	RMSHE	RSIHE	ADPHE	Ours
‘residence’	1.0765	1.0791	1.0795	1.0529	1.0578	1.0664	1.0423	1.0409
‘car’	0.9312	0.9116	0.9393	0.9045	0.9044	0.8591	0.8935	0.8580
‘guardrail’	0.8900	0.8888	0.8965	0.9043	0.9054	0.9222	0.9354	0.8687
‘building’	1.0108	1.0083	1.0100	1.0059	1.0092	0.9949	0.9894	0.9388
‘chair’	1.0219	1.0143	1.0232	1.0136	1.0089	1.0188	1.0093	1.0083
‘table’	0.5149	0.5462	0.5322	0.6583	0.6733	0.6725	0.7245	0.5060
‘leg’	0.9830	0.9247	0.9860	0.9771	0.9067	0.9219	0.8301	0.6563
‘computer’	0.7749	0.7911	0.7917	0.8895	0.8776	0.8826	0.9063	0.7720
Ave.	0.9004	0.8955	0.9073	0.9257	0.9179	0.9173	0.9163	0.8311

,

Table 7. Quantitative comparisons of DE for test images, and bold values indicate the top two results.

	GHE	BBHE	MMBEBHE	DISHE	RMSHE	RSIHE	ADPHE	Ours
‘residence’	4.6137	4.7407	4.7440	4.6254	4.7362	4.5887	4.7578	4.7814
‘car’	5.3755	5.6892	5.6620	5.6261	5.6718	5.6029	5.3767	5.7270
‘guardrail’	5.5252	5.7597	5.7517	5.6302	5.6942	5.7341	5.4055	6.0788
‘building’	5.0128	5.2030	5.2001	5.1155	5.1643	5.0980	5.1544	5.5099
‘chair’	5.4534	5.6513	5.6509	5.5331	5.5749	5.4956	5.4883	5.6769
‘table’	5.1794	5.3335	5.3334	5.2613	5.3089	5.2492	5.1603	5.6247
‘leg’	4.9969	5.4267	5.3248	5.2583	5.4086	5.2907	5.3053	5.4650
‘computer’	4.0868	4.1793	4.1791	3.9492	4.1579	4.0408	4.1892	4.4657
Ave.	5.0304	5.2479	5.2307	5.1248	5.2146	5.1375	5.1046	5.4161

and

Table 8. Quantitative comparisons of AME for test images, and bold values indicate the top two results.

	GHE	BBHE	MMBEBHE	DISHE	RMSHE	RSIHE	ADPHE	Ours
‘residence’	5.8276	6.0680	6.5956	13.6031	12.9730	16.2548	16.0403	13.7693
‘car’	6.7853	6.6191	7.4672	14.3339	14.2674	15.9727	15.3144	15.4704
‘guardrail’	4.0183	4.1726	4.5488	10.1015	9.9692	11.4671	10.8266	11.5326
‘building’	5.3333	5.4193	6.0829	12.6270	12.4008	14.1731	13.8817	14.4274
‘chair’	9.1302	8.9374	9.7847	17.0440	16.8935	18.8566	17.3644	17.4060
‘table’	8.8114	10.1058	9.5605	15.6672	15.5742	17.0279	17.8687	19.1004
‘leg’	7.0446	6.1563	7.7684	15.353	16.7957	17.6740	15.1999	16.9298
‘computer’	4.4752	5.1033	5.4196	11.2161	10.5654	13.6330	17.8573	17.8747
Ave.	6.4282	6.5727	7.1534	13.7432	13.6799	15.6324	15.5441	15.8138

list the statistical results of the six indexes and the average values are also given in the last row of each table. To present a clearer comparison, the corresponding line graphs are also presented in Figure 20.

On the one hand, $\mathsf{\eta }$ denotes the target/background contrast, and it also measures the quantity of clutters to a certain degree. As is shown in Table 3 and Figure 20a, our method achieves the best average $\mathsf{\eta }$, implying that the foregrounds become brighter while the background clutters become darker after our enhancement. More specifically, we argue that it is the strategy of sub-histogram segmentation that enlarges the intensity difference between the fore- and background. We notice that BBHE, GHE, and MMBEBHE also result in low $\mathsf{\eta }$ values (<2) in most scenes, which matches the experimental fact that the intensity differences between target and background regions are always remarkable in their enhanced images. Besides, the average $\mathsf{\eta }$ of RSIHE is the worst, which is about 10% larger than ours and it proves that the serious over-enhancement is somewhat related to the degradation of global contrast.

Unlike its poor performance related to $\mathsf{\eta }$, RSIHE has the most satisfactory PSNR on the whole, while our method still takes second place in this metric. Since PSNR reflects the degree of image degradation, we consider that if the PSNR value is excessively large, a large amount of detail information is lost at the same time. That is the reason why the enhanced images processed by RSIHE are fogged and their PSNRs are always very high. By contrast, our method maintains a balance between PSNR (slightly larger than ADPHE and RMSHE) and exhibits the preservation of details, meaning that not only the visual performance presented in Section 4.2, but also the quantitative evaluation of image degradation, is qualified.

As is mentioned above, $\mathsf{\rho }$ is a comprehensive index measuring the overall performance of image enhancement. Since our method results in comparatively remarkable performances in both $\mathsf{\eta }$ and PSNR, its $\mathsf{\rho }$ values of all the test images thus keep the advantage. Relying on its outstanding achievements in $\mathsf{\eta }$, the average $\mathsf{\rho }$ value of our method is almost two times that of RSIHE, which possesses the second largest value. One point that should be noticed is that those algorithms overly emphasizing input brightness preservation, e.g., BBHE, MMBEBHE, DISHE, RMSHE and RSIHE, sacrifice the fore- and background contrast instead, which leads to their poor performances in $\mathsf{\eta }$ and $\mathsf{\rho }$.

$\mathsf{\Omega }$ is commonly utilized to reflect the smoothness of image. As is shown in Figure 20d, BBHE and our method obtain the top two $\mathsf{\Omega }$ values, whereas DISHE presents the worst one. Since the background occupies the majority of an image and the background histograms are not equalized in our method, the noise is not amplified and the smoothness is thus satisfactory. However, some grayscale histogram equalization-based algorithms, like GHE and MMBEBHE, quite easily generate over-enhancement in high grayscales and image noise precisely tends to center on high grayscales, so the effects of the noise suppression of these algorithms greatly decrease. According to Table 6, we also see that ADPHE’s capability of removing noise through the adaptive upper threshold is limited in IR images considering its comparatively high $\mathsf{\Omega }$ values.

DE is a significant index that directly measures the degree of over-enhancement. Through observing its definition, we can easily find that if more grayscales are merged by equalization, the entropy of the image will decrease more. Clearly, our method has a distinct advantage in DE which is approximately 3% greater than the second largest one (BBHE). On the contrary, GHE and DISHE decrease the DE value of the original image seriously, and this is mainly because these methods tend to merge different gray levels together through enhancement. Interestingly, RSIHE achieves a satisfactory performance in PSNR, but its DE performance is poor. Furthermore, BBHE undergoes an opposing situation. We infer that a high DE value can indicate relatively low image degradation to a certain degree.

To evaluate the local contrast of the enhanced image, Table 8 reports the comparison result in terms of AME. Among all the comparative algorithms, the average AME of our method is the best, while the one of GHE is the worst. Even though the $\mathsf{\eta }$ values of GHE and MMBEBHE, which represent the global contrast, are large, plenty of local regions turn to be homogenous due to the fact that it is highly likely that pixels in the same region will merge in these algorithms, contributing to the low AMEs.

4.3.2. Evaluation of Running Time

A comparison of running time is made in

Table 9. Comparison of execution time with conventional methods (s), and bold values indicate the top two results.

	GHE	BBHE	MMBEBHE	DISHE	RMSHE	RSIHE	ADPHE	Ours
residence	0.1287	0.0112	0.1253	0.0165	0.1336	0.1396	0.0206	9.5613
car	0.1272	0.0106	0.1223	0.0161	0.1387	0.1405	0.0269	9.7212
guardrail	0.1255	0.0123	0.1225	0.0159	0.1337	0.1415	0.0207	9.4910
building	0.1282	0.0105	0.1269	0.0157	0.1367	0.1474	0.0225	9.3013
room	0.1288	0.0122	0.1264	0.0151	0.1356	0.1445	0.0216	8.9021
cabinet	0.1278	0.0113	0.1235	0.0154	0.1338	0.1493	0.0246	9.3456
leg	0.1283	0.0119	0.1247	0.0165	0.1335	0.1409	0.0248	9.5722
computer	0.1297	0.0110	0.1267	0.0174	0.1377	0.1494	0.0277	8.9022
Ave.	0.1280	0.0114	0.1248	0.0161	0.1354	0.1441	0.0237	9.3496

. All the algorithms are implemented for five times and the average running time is recorded.

As we can see, BBHE achieves the fastest running speed, which is about 88 fps. Besides, DISHE and ADPHE also perform well, and their running speeds can reach more than 60 fps and 40 fps, respectively. In regard to our method, the PSO-based optimization consists of multiple iterations of the procedures introduced in Section 3.2. That is to say, these parts of codes need to be executed for $\mathrm{N}\times {\mathrm{t}}_{\max }$ times in total, which leads to the long computation time. As a matter of fact, computation burden is a common problem existing in all the iterative algorithms which needs to be urgently settled. Motivated by the extensive applications of multi-core digital signal processors, e.g., GPU, in practical engineering, we propose to increase our running efficiency using the strategy of parallel processing. Considering that each particle in PSO is individual, the executable codes of all the particles in the same iteration can be regarded as a group of parallel tasks and can be implemented in different cores at the same time. Under the circumstances, even real-time running can be expected.

In order to demonstrate the advantage of PSO optimization among concurrent methods, three other concurrent optimization algorithms: the genetic algorithm (GA) [44], ant colony optimization (ACO) [45], and the bat algorithm (BA) [46] are applied to make a comparison. For a fair comparison, all these algorithms are implemented with the same number of particles and iterations ($\mathrm{N}=10,{\mathrm{t}}_{\max }=10$), and

Table 10. Comparison of running time with concurrent algorithms (s), and bold values indicate the top results.

	GA	ACO	BA	PSO
residence	10.5096	10.8130	9.3390	9.5613
car	10.9597	10.5610	9.5279	9.7212
guardrail	10.7847	10.5709	9.5572	9.4910
building	10.0764	10.5479	9.3078	9.3013
room	10.8000	10.8842	9.0858	8.9021
cabinet	10.7875	10.2348	9.4933	9.3456
leg	10.9137	10.8163	9.4611	9.5722
computer	10.5576	10.1879	8.9899	8.9022
Ave.	10.6737	10.5770	9.3453	9.3496

clearly reports the running time of the four algorithms. Note that the execution time listed is the average of five times of running.

As is indicated from Table 10, the average running time of GA and ACO is obviously longer than that of BA and PSO. This is because GA contains several genetic operations, e.g., crossover and mutation, and ACO contains a series of combinatorial optimizations. BA and PSO achieve almost the same average execution time due to the fact that their update mechanisms of position and velocity are similar.

4.3.3. Evaluation of Mean Brightness Error Ratio

Intended to further verify the advantage of our method when it is compared with other concurrent methods, a comparison of the mean brightness error ratio $\mathsf{\delta }$, which reflects the optimization ability, is made in this section. Here, the definition of $\mathsf{\delta }$ is given below:

$$\mathsf{\delta }=\frac{\left|\mathrm{mean}(\mathrm{O})-\mathrm{mean}(\mathrm{Ref})\right|}{\mathrm{mean}(\mathrm{Ref})},$$

(43)

where, mean(O) and mean(Ref) stand for the mean intensity values of the enhanced image and the reference image, respectively.

Table 11. Optimization results of ${\mathrm{m}}_1^{\prime }$.

	Residence	Car	Guardrail	Building	Room	Cabinet	Leg	Computer
GA	60	38	21	52	4	32	55	0
ACO	55	41	22	55	3	32	56	0
BA	61	38	24	57	9	28	65	0
PSO	56	41	21	54	3	31	54	0

records the optimized values of ${\mathrm{m}}_1^{\prime }$ and

Table 12. Mean brightness error ratio $\mathsf{\delta }$, and bold values indicate the top results.

	Residence	Car	Guardrail	Building	Room	Cabinet	Leg	Computer	Ave.
GA	0.0276	0.0209	0.0016	0.0158	0.0069	0.0108	0.0099	0.1435	0.0296
ACO	0.0121	0.0077	0.0058	0.0081	0.0031	0.0108	0.0147	0.1435	0.0257
BA	0.0355	0.0209	0.0207	0.0219	0.0400	0.0182	0.0854	0.1435	0.0483
PSO	0.0041	0.0077	0.0016	0.0019	0.0031	0.0019	0.0012	0.1435	0.0206

records the mean brightness errors. It can be observed from these two tables that all four of the concurrent algorithms are able to generate reasonable optimization results and the average values of $\mathsf{\delta }$ are all lower than 3.5. However, we need to notice that apart from the ‘computer’ image, the $\mathsf{\delta }$ values of our method are superior than most of the others, which convincingly verifies the optimization ability of our algorithm. We consider that this is mainly owing to the memory mechanism of PSO. Furthermore, the reason why the $\mathsf{\delta }$ values of the ‘computer’ image are relatively large is that the brightness of the input itself is much higher than the reference image, so the mean brightness of the enhanced image cannot approach the reference even if ${\mathrm{m}}_1^{\prime }=0$.

5. Conclusions

In this paper, a new IR image enhancement method using adaptive histogram segmentation and local contrast weighted distribution is presented. In order to separately address the fore- and background, the raw grayscale histogram is firstly smoothed by Gaussian filtering and LOWESS and partitioned via local minima. Then, a metric called grayscale density is developed to recognize whether each sub-histogram belongs to the foreground or background. For the foreground sub-histograms, a local contrast weighted distribution is calculated via modified sigmoid transformation and is applied to replace the intensity distribution, after which the detail information can be enhanced by equalization based on the proposed distribution. For the background sub-histograms, a strategy where the corresponding proportions of the whole dynamic range are maintained is considered to avoid over-enhancement. Also, an exponential factor is designed based on a characteristic of human vision and is employed to adjust the re-mapping ranges for the purpose of decreasing the effect of noise. Finally, the overall brightness of the output is corrected through the PSO algorithm based on a reference image with an appropriate mean brightness. A number of both qualitative and quantitative experiments prove that our method outperforms other conventional methods in different scenes.

We would like to acknowledge that there are 13 parameters utilized in our algorithm and most of the values are empirically set (see Section 4.1). Fairly speaking, this number of parameters may be more than other existing methods, although it can achieve a more satisfactory performance in IR image enhancement. In future work, we will pay more attention to this issue and improve the robustness and simplicity of our algorithm further.

Also, we will concentrate on transplanting the presented algorithm to an appropriate hardware platform. Compared with the state-of-the-art methods, the PSO iteration in our method hugely increases the computational complexity. Under the circumstances that we would like to improve the accuracy of optimization, more particles and more iterative rounds are needed, which will further increase the running time. However, we consider that this problem can be solved if multi-core processors, e.g., GPU, are utilized, because every particle can be regarded as an individual task and be implemented in different cores simultaneously.