Histogram Equalization Algorithm Based on Sobel Gradient and Its Application on Infrared Images
-
摘要: 为了能在动态范围压缩的同时增强红外图像的对比度，提出了一种基于Sobel梯度直方图均衡算法（gradient histogram equalization，GHE）。与以往的直方图均衡化（histogram equalization，HE）方法不同，该方法自适应地为图像强梯度的灰阶分配高对比度，保留并增强16 bit图像中更多的细节。随后使用双Gamma映射对映射曲线进行调整，有效地抑制图像亮部的过曝现象，同时提高暗部的细节。该方法相比于传统的直方图均衡化算法在暗区细节处理、过曝抑制、对比度增强等方面都有较好的效果。Abstract: A Sobel gradient histogram equalization(GHE) algorithm is proposed to enhance the contrast of infrared images during dynamic range compression. In contrast to previous histogram equalization(HE) methods, this method adaptively assigns a high contrast to the strongly graded parts of the image, preserving and enhancing more details in the 16-bit image. Dual Gamma mapping is then used to adjust the mapping curve to effectively suppress overexposure in the bright parts of the image while improving the detail in the shadows. Compared with the traditional histogram equalization algorithm, this method has better effects on dark area detail processing, overexposure suppression, and contrast enhancement.
-
Keywords:
- infrared image /
- gradient histogram equalization /
- dual Gamma mapping
-
0. 引言
红外图像处理技术广泛应用于军事、医疗、防灾等领域。其主要应用是提高红外成像设备的成像质量，处理现有的红外图像。红外图像的成像特征来自目标的能量辐射，其特点是能量不均匀，分辨率对比度较低，亮度分布不均[1]。因此，红外图像对比度增强的方法在实际应用中具有重要意义[2]。
图像的对比度增强是图像信号处理中流行的一个研究方向。其主要目的是扩大图像对比度，提高图像的视觉效果，反映一些隐藏在图像中的细节[3]。目前常见的对比度增强方法是针对可见光图像进行增强，大致分为两类：①直接增强输入图像[4]。②通过间接的方法增强输入图像[5]。通常直接增强很容易对噪声进行过度的增强。间接增强法则利用输入灰度直方图的统计量，对灰度直方图进行处理，最后通过处理后的直方图建立输入图像的灰度与输出图像的灰度之间的映射关系，以增强图像的对比度。
基于直方图均衡化的图像对比度增强方法，在可见光图像处理领域得到了广泛的应用，在红外图像对比度增强领域也具有较好的效果。对于给定的图像，该方法使用图像的累积分布函数来扩展图像的动态范围，以提高图像的对比度[6]。然而不同灰度的增强程度与其直方图的分布有关，因此直方图均衡化的方法容易产生灰度平均位移的问题，导致图像目标或噪声被严重衰减或者过度增强[7]。多位学者对基于直方图均衡化的方法进行了研究。Hum Y. C.等人使用全局加局部的对比度增强算法，增强图像的低对比度区域，抑制高对比度区域[8]。J. Bae等人先根据对数域的直方图将图像分成K类，然后每一个类使用不同的log域映射算法[9]。Weitao Deng等人通过递归将图像的直方图分成多个子直方图，通过特定的方式加权平均得到均分分布的直方图，最后根据原始的累计直方图将两种直方图加权得到最终用于直方图均衡的直方图[10]。T. Arici等人提出多种直方图修改方案，根据其特征提出了一种新的基于局部对比度的直方图统计方法，在小计算量下达到较好的效果[11]。
本文提出一种基于Sobel梯度的直方图均衡算法，可以获得更好的图像增强效果，同时通过双Gamma映射调整映射曲线，可有效避免梯度过大导致的亮区过曝、暗区欠曝的现象。
1. 提出的算法
1.1 直方图均衡化
红外探测器的输出数据一般为16 bit的高动态范围图像，通过对16 bit图像进行直方图均衡化，可以得到8 bit的灰度图像，同时实现了对红外图像的动态范围压缩和增强[12]。假设红外探测器输出的某一帧16 bit图像共有L(1≤L≤216)个灰度等级，则图像的直方图可定义如下：
$$\operatorname{hist}(l)=n(l), \quad l_{\min } \leqslant l \leqslant l_{\max }$$ (1) 式中：n(l)表示灰阶为l的像素点在全图出现的次数；lmin为全图最小灰阶；lmax为全图最大灰阶。图像直方图描述了所有灰度级在全图范围内的分布情况，直方图均衡化的核心思想是通过一种特定的映射关系将原图直方图分布中较为集中的像素映射成均匀分布[13]。
将直方图正则化后，可得到图像的离散概率密度函数，定义如下：
$$ p\left( l \right) = \frac{{{\rm{hist}}\left( l \right)}}{N} $$ (2) 式中：p(l)表示灰阶为l的像素点在全图出现的概率；N为全图像素总个数。由此，直方图均衡化过程可描述如下：
$$p^{\prime}=T(p)$$ (3) 式中：p′表示直方图均衡化后的离散概率密度函数，为均匀分布，即所有灰阶出现的概率相同。通过求解上式，可得到映射关系T为累计分布函数，可定义如下：
$$ c\left( l \right) = \sum\limits_{i = {l_{\min }}}^l {p\left( i \right)} $$ (4) 式中：c(l)∈[0, 1]，将其乘以8 bit最大灰度值255后即可得到16 bit图像到8 bit图像的映射关系。对于16 bit图像上某一个灰度值为l的像素，映射后其在8 bit图像上的对比度可描述如下：
$$ \Delta l = 255 \times \left( {c\left( l \right) - c\left( {l - 1} \right)} \right) = 255 \times p\left( l \right) $$ (5) 从上式可以看出，某一个灰阶在8 bit图像的亮度及对比度与其在16 bit上出现的概率有关。其出现概率越高，映射后与相邻灰度级所占用的8 bit动态范围越大，表现为图像的对比度越大。当某个范围内的灰阶在全图所占比例很大时，映射后会占用8 bit图像很大的动态范围，从而压缩其余部分所能占用的动态范围。但占比很大的像素并不一定会存在很多细节，从而导致占比少、细节多的部分图像细节被严重压缩。因此本文提出一种基于Sobel梯度的直方图均衡算法，可以自适应地为图像强梯度的灰阶分配高对比度，能有效地保留并增强原图的细节。针对梯度累计过大导致的梯度直方图峰值过高，从而出现的对比度分配过度现象，使用一种双Gamma映射对最终的映射曲线进行调整，可有效地抑制图像亮部的过曝现象，同时提高暗部的细节。
1.2 梯度直方图均衡化
对于输入的16 bit图像srcImage，首先计算其Sobel梯度图像gradImage，梯度图像中每一个点计算公式如下：
$$ \begin{array}{l} gx\left( {i, j} \right) = \left( {{\rm{srcImage}}\left( {i + 1, j - 1} \right) + 2 \times {\rm{srcImage}}\left( {i + 1, j} \right) + {\rm{srcImage}}\left( {i + 1, j + 1} \right)} \right) - \hfill \\ \quad \quad \quad \, \, \, \, \,\,\,\, \left( {{\rm{srcImage}}\left( {i - 1, j - 1} \right) + 2 \times {\rm{srcImage}}\left( {i - 1, j} \right) + {\rm{srcImage}}\left( {i - 1, j + 1} \right)} \right) \hfill \end{array} $$ (6) $$ \begin{array}{l} gy\left( {i, j} \right) = \left( {{\rm{srcImage}}\left( {i + 1, j + 1} \right) + 2 \times {\rm{srcImage}}\left( {i, j + 1} \right) + {\rm{srcImage}}\left( {i - 1, j + 1} \right)} \right) - \hfill \\ \quad \quad \quad \, \, \, \, \,\,\,\, \left( {{\rm{srcImage}}\left( {i + 1, j - 1} \right) + 2 \times {\rm{srcImage}}\left( {i, j - 1} \right) + {\rm{srcImage}}\left( {i - 1, j - 1} \right)} \right) \hfill \end{array} $$ (7) $$ {\rm{gradImage}}\left( {i, j} \right) = \sqrt {gx{{\left( {i, j} \right)}^2} + gy{{\left( {i, j} \right)}^2}} $$ (8) 式中：(i, j)为图像坐标，对梯度图像进行归一化：
$$ {\rm{gradImageNor}}\left( {i, j} \right) = {\rm{gradImage}}\left( {i, j} \right)/\max \left( {{\rm{gradImage}}} \right) $$ (9) 式中：max(gradImage)为梯度图像最大值。
红外图像在经过非均匀性校正后依然会存在一定程度的非均匀性噪声，特别是在局部动态范围很低的背景上表现的尤为明显。非均匀性噪声往往会产生很小的梯度值，但当噪声数量较多时，对小梯度进行累计，依然会在梯度直方图上产生较高的bin值，因此需要在梯度统计时进行筛选，避免为非均匀性噪声产生较高的bin值。本文通过设置固定的阈值gradTh对梯度进行筛选，将小于gradTh的梯度值赋值为0，防止后续统计直方图时对其进行累计。
$$ \begin{array}{l} {\rm{gradImageTh}}\left( {i, j} \right) = \hfill \\ \left\{ {\begin{array}{*{20}{c}} 0&{{\rm{gradImageNor}}\left( {i, j} \right) < {\rm{gradTh}}} \\ {{\rm{gradImageNor}}\left( {i, j} \right)}&{{\rm{gradImageNor}}\left( {i, j} \right) \geqslant {\rm{gradTh}}} \end{array}} \right. \hfill \end{array} $$ (10) 在计算梯度图像后，需要计算输入图像的梯度直方图，计算过程与正常直方图统计类似，只是将灰度级出现的个数累计更改为对应位置的梯度累计，16 bit图像中某个灰度级k的梯度直方图bin值histGrad(k)计算公式如下：
$$ {\rm{histGrad}}\left( k \right) = \sum\limits_{i = 1}^h {\sum\limits_{j = 1}^w {f\left( {i, j, k} \right) \times g\left( {i, j} \right)} } $$ (11) 式中：h、w分别为图像高和宽；k表示16 bit图像中出现的灰度级，lmin≤k≤lmax。当srcImage(i, j)＝k时，f(i, j, k)＝1，否则f(i, j, k)＝0。对梯度直方图累计，得到其累计直方图，计算公式如下：
$$ {\rm{histCum}}\left( k \right) = \sum\limits_{l = {l_{\min }}}^k {{\rm{histGrad}}\left( l \right)} $$ (12) 对累计直方图归一化即可得到最终的映射曲线：
$$ {\rm{map}}\left( i \right) = \frac{{{\rm{histCum}}\left( i \right)}}{{\max \left( {{\rm{histCum}}} \right)}} $$ (13) 式中：lmin≤k≤lmax，max(histCum)表示累计直方图最大值，即histCum(lmax)。同理，将其乘以8 bit最大灰度值255后即可得到16 bit图像到8 bit图像的映射关系。
将普通直方图和梯度直方图均归一化后，在某一场景下对比如图 1所示，其中(a)为原始图像归一化到0~1后直接显示，(b)为该图像的直方图和梯度直方图。(b)中红色虚线框区域表示天空所处灰度级范围，绿色虚线框区域表示树木所处灰度级范围。在原图中，由于天空和树木占比很大，所以直方图中在天空所处的灰度级和树木所处的灰度级范围产生了很高的bin（指直方图中横坐标的某一个范围，其高度代表图像灰阶中处于该范围的像素点个数，本文中横坐标范围固定为1个灰阶），其余部分的bin很低。而在梯度直方图中，由于天空部分没有很高的梯度值，所以bin被拉低，而树木存在较强的梯度值高bin得到了保留，图像中建筑等其余部分由于存在很高的梯度，bin得到了拉伸。
1.3 双Gamma映射
红外图像的观察者通常更加关注温度明显异于环境温度的目标，如导弹尾焰、猎物、房屋漏水处等，这些目标在16 bit数据中的灰阶分布一般都处于全图灰阶的最大值或最小值附近，即对应为图像的亮区和暗区。梯度直方图在进行梯度累计时，直方图bin的高度不仅与梯度值相关，还与bin对应的灰度级在图像中的个数有关。图像背景一般会占到图像绝大部分，导致背景灰度等级对应的bin会产生较大的值，从而导致映射后背景会占用很大的8 bit动态范围。而在图像的亮区和暗区即使存在很大的梯度，由于像素个数较少，其梯度直方图的bin依然不会很高，且8 bit动态范围被背景占用了很大一部分，导致映射后的图像亮区和暗区可以使用的动态范围很小，从而出现亮区过曝、暗区欠曝的现象。如图 2所示，场景中人物头部基本都发生了过曝。
针对上述问题，本文引入一种双Gamma映射曲线[14]对最终得到的map进行二次映射。修正后的映射曲线可计算如下：
$$ {\rm{ma}}{{\rm{p}}_{{\rm{cor}}}}\left( i \right) = \frac{{{\rm{map}}{{\left( i \right)}^{1/{\rm{gamma}}}} + \left( {1 - {{\left( {1 - {\rm{map}}\left( i \right)} \right)}^{1/{\rm{gamma}}}}} \right)}}{2} $$ (14) 式中：gamma为手动设置的值，可根据场景手动调整，不同gamma值下的反S型曲线如图 3所示。当gamma＝1时，相当于对原始曲线不做映射。随着gamma值的增大曲线两端的斜率会增加，曲线中部的斜率会降低，意味着图像亮区和暗区的对比度被拉伸，图像灰度级中间区域对比度被压缩。
2. 实验与结果分析
为验证本文提出算法的有效性，在非制冷1280×1024分辨率探测器和制冷640×512分辨率探测器采集的raw数据上进行了仿真验证。其中非制冷1280×1024分辨率探测器的raw数据在进行动态范围压缩前进行了非均匀性校正和去横竖纹处理，制冷640×512分辨率探测器的raw数据进行了非均匀性校正和去噪处理。
2.1 梯度直方图均衡化效果分析
如图 4所示，对比了室外场景下直方图均衡化和梯度直方图均衡化的效果，以及两种算法获得的映射曲线。图(a)为直方图均衡化效果，图(b)为梯度直方图均衡化效果，图(c)为两种算法获得的映射曲线。其中图(c)中红色虚线框表示直方图均衡化映射曲线对应图像中红色区域天空部分，绿色虚线框表示直方图均衡化映射曲线对应图像中绿色区域建筑部分。可以看到在该场景中，由于天空占比比较大，因此原本在16 bit图像中动态范围很小的天空部分，在经过直方图均衡化算法后，在8 bit图像中占用了很大的动态范围内。从而导致天空部分对比度被过度拉伸，增强了天空部分的噪声。另一方面，由于天空部分占用了很多动态范围，原本在16 bit图像中动态范围很大的建筑部分可以被分配的动态范围变得很少，导致建筑部分对比度很低，部分细节无法被凸显。由于梯度直方图根据像素梯度进行8 bit图像动态范围的分配，天空部分梯度值很小，从而在8 bit图像中被分配了很少的动态范围，抑制了噪声；建筑部分梯度值大，从而在8 bit图像中被分配了很大的动态范围，凸显了建筑细节。
为量化分析对比度的分配情况，本文定义对比度计算公式如下：
$$ c = \sum\limits_\delta ^{} {{\delta ^2} \times {P_\delta }} $$ (15) 式中：δ表示图像中任意相邻两点的灰度差，本文计算时取四领域，即只计算水平方向和垂直方向相邻的像素；Pδ表示灰度差为δ在整个统计过程中出现的概率。如表 1所示，计算了图 4中红色框区域天空部分和绿色框区域建筑部分的局部对比度以及全图对比度。可以看出梯度直方图均衡化将非显著性区域的天空部分对比度从原有的3.74压制到了0.19，而为显著性区域的建筑部分分配了更高的对比度，从之前的54.33提升到了123.96，同时全图的对比度也有明显增强。
表 1 两种算法对比度对比Table 1. Comparison of HE with GHEAlgorithms Sky Architecture Picture HE 3.74 54.33 74.14 GHE 0.19 123.96 98.19 2.2 噪声抑制效果分析
如图 5所示，对比了室内某一场景下的不进行噪声抑制和增加噪声抑制的效果及映射曲线，其中阈值gradTh＝0.05。图(a)为直方图均衡化效果，图(b)为梯度直方图均衡化效果，图(c)为噪声抑制的梯度直方图均衡化效果，图(d)为3种算法获得的映射曲线。其中图(d)中红色虚线框表示直方图均衡化映射曲线对应图像中背景区域。可以看到在该场景下，原本的直方图均衡化和梯度直方图均衡化均会放大背景的噪声。其主要原因为背景噪声的梯度值虽然很小，但数量很多，经过累计后背景所在的灰度级依然会在直方图均衡上产生很高的bin，导致其在8 bit图像中占用很大的动态范围，从而增强噪声。但增加阈值限制后，背景所在灰度级梯度值不会被累计，从而避免该问题，同时增强目标的对比度。但同时与背景噪声灰度级接近的细节也会被抑制，需要根据场景适当调整该参数，图中gradTh＝0.05。
2.3 双Gamma映射效果分析
如图 6和图 7所示，分别对比了室外和室内场景下的不进行双Gamma映射和增加双Gamma映射的效果及映射曲线，其中gamma＝2.2。图(a)为直方图均衡化效果，图(b)为梯度直方图均衡化效果，图(c)为双Gamma映射的梯度直方图均衡化效果，图(d)为3种算法获得的映射曲线。图 6中(a)~(c)红色框选区域与(d)红色虚线框选的直方图均衡化映射曲线相对应。可以看到在该场景下，直方图均衡化和梯度直方图均衡化均在红色框选区域内产生了很严重的过曝，主要由于原图中的树木和远处建筑所处灰度等级占全图像素的绝大部分，且该部分灰阶也具备一定的梯度值，进行梯度累计后产生很高的bin，导致直方图均衡和梯度直方图均衡都会为该部分在8 bit图像上分配很大的动态范围，从而压缩了其他部分可用的动态范围，导致反射太阳的楼体出现过曝，看不到窗户的细节。增加双Gamma映射后，16 bit图像动态范围的中间部分映射后的动态范围被压缩，过曝区域映射后的动态范围被拉伸，使得窗户的细节被保留。
同样的，图 7中(a)~(c)红色框选区域与(d)红色虚线框选的双Gamma映射的梯度直方图均衡化映射曲线相对应。该场景下直方图均衡化和梯度直方图均衡化效果，人物头部都发生了过曝，细节损失严重，增加双Gamma映射后，人物头部的映射曲线会更加陡峭，意味着分配了更高的对比度，从而使得人物头部细节凸显。
3. 结语
本文提出一种基于Sobel梯度的直方图均衡化算法，可以根据全图的梯度分布情况，自适应地为强梯度的灰阶分配高对比度。同时使用一种双Gamma映射对最终的映射曲线进行调整，使得红外图像动态范围压缩过程中，不会出现明显的过曝和欠曝现象，能更好地保留亮区和暗区的细节。
但本文提出的算法仍有一定的局限性。算法在自适应为强梯度的灰阶分配高对比度的同时，不可避免会导致图像弱纹理的对比度被压缩，从而导致一定的细节丢失。另外，在低动态范围场景，动态范围压缩的过程逆转为动态范围增强，使用本文算法较容易出现图像失真现象。
在后续的研究中，可以尝试将图像的梯度信息和其他信息（场景动态范围、局部空间熵等）进行融合，从多个维度考虑对比度的自适应分配问题，避免出现图像失真。同时，也可以考虑结合局部对比度增强的方法，避免弱纹理被压缩的问题。
-
表 1 两种算法对比度对比
Table 1 Comparison of HE with GHE
Algorithms Sky Architecture Picture HE 3.74 54.33 74.14 GHE 0.19 123.96 98.19 
下载: 导出CSV
-
[1] LIU C, SUI X, KUANG X, et al. Optimized contrast enhancement for infrared images based on global and local histogram specification[J]. Remote Sensing, 2019, 11(7): 849. DOI: 10.3390/rs11070849
[2] 李凌杰, 陈菲菲. 基于改进直方图的红外图像增强方法[J]. 航空兵器, 2022, 29(2): 101-105. https://www.cnki.com.cn/Article/CJFDTOTAL-HKBQ202202015.htm LI Lingjie, CHEN Feifei. Infrared image enhancement method based on improved histogram[J]. Aero Weaponary, 2022, 29(2): 101-105. https://www.cnki.com.cn/Article/CJFDTOTAL-HKBQ202202015.htm
[3] CHEN X, LV L. A compositive contrast enhancement algorithm of IR image[C]//International Conference on Information Technology & Applications of IEEE, 2013, DOI: 10.1109/ITA.2013.20
[4] LIU N, CHEN X. Infrared image detail enhancement approach based on improved joint bilateral filter[J]. Infrared Physics & Technology, 2016, 77: 405-413.
[5] HUANG Z, ZHANG T, LI Q, et al. Adaptive gamma correction based on cumulative histogram for enhancing near-infrared images[J]. Infrared Physics & Technology, 2016, 79: 205-215.
[6] JUNG C, LU W. Tone-preserving contrast enhancement in images using rational tone mapping and constrained optimization[C/OL]. 2016 IEEE Visual Communications and Image Processing, https://ieeexplore.ieee.org/document/7805525.
[7] LI S, JIN W, LI L, et al. An improved contrast enhancement algorithm for infrared images based on adaptive double plateaus histogram equalization[J]. Infrared Physics & Technology, 2018, 90: 164-174.
[8] Hum Y C, Tee Y K, W. -S. Y, et al. A contrast enhancement framework under uncontrolled environments based on just noticeable difference[J]. Signal Processing, 2022, 103: 116657
[9] Bae J, Kim K, Yun Y J, et al. Adaptive tone-mapping operator for HDR images based on image statistics[C]//Proc. of IEEE, 2011, DOI: 10.1109/TENCON.2011.6129047
[10] DENG Weitao, LIU Lei, CHEN Huateng, et al. Infrared image contrast enhancement using adaptive histogram correction framework[J]. Optik, 2022, 271: 170114. DOI: 10.1016/j.ijleo.2022.170114
[11] Arici T, Dikbas S, Altunbasak Y. A histogram modification framework and its application for image contrast enhancement[J]. IEEE Trans Image Process, 2009, 18(9): 1921-1935. DOI: 10.1109/TIP.2009.2021548
[12] Ok J, Lee C. HDR tone mapping algorithm based on difference compression with adaptive reference values[J]. Journal of Visual Communication and Image Representation, 2016, 43: 61-76.
[13] Arici T, Dikbas S, Altunbasak Y. A histogram modification framework and its application for image contrast enhancement[J]. IEEE Trans Image Process, 2009, 18(9): 1921-1935. DOI: 10.1109/TIP.2009.2021548
[14] Lee S. An efficient content-based image enhancement in the compressed domain using Retinex theory[J]. IEEE Transactions on Circuits & Systems for Video Technology, 2007, 17(2): 199-213.
-
期刊类型引用(0)
其他类型引用(1)
计量
- 文章访问数: 164
- HTML全文浏览量: 40
- PDF下载量: 65
- 被引次数: 1
