计算机视觉入门之图像处理＜四＞：图像直方图均衡化

往期文章回顾：
计算机视觉入门之<零>
计算机视觉入门之图像处理<一>：图像处理基础概念
计算机视觉入门之图像处理<二>：图像处理基础概念
计算机视觉入门之图像处理<三>：图像插值方法

---

前言

掌握了图像处理基础知识，该系列文章进一步了解图像处理基本方法，本篇文章主要包括图像直方图理论、直方图均衡化原理以及代码实现等方面。

灰度直方图基本概念

灰度范围为$\left[0,L-1\right]$的数字图像中，将其中每一个灰度级看成一个随机变量，则随机变量的取值情况就表示了该图像的统计特性，改用直方图来描述该特性，则称之为灰度直方图。
灰度直方图是灰度级的离散化函数，表示为：
$$h(r_k)=n_k$$其中$r_k$表示第$k$级灰度值，$n_k$表示像素为$r_k$的像素个数。
灰度直方图的归一化：$$p(r_k)=\frac{n_k}{MN}$$式中$MN$表示图像的行数和列数的乘积，归一化后表示的意义是灰度值$r_k$在图像中出现概率的估计。归一化后的直方图所有分量满足下列关系式：$$\sum _{k=0}^{L-1}p(r_k)=1$$

如图所示的四个基本灰度级的特征花粉图像，分别表现为暗、亮、低对比度和高对比度，灰度图横轴表示灰度值$r_k$， 纵轴表示$n_k$或归一化后的值$p(r_k)$。由图可以得出，直方图越靠近灰度值小的一端，图像越暗，反之；直方图集中在灰度级中部的图像具有较低的对比度，而高对比度的图像直方图几乎覆盖了所有灰度级，而且各个灰度值的分布较为均匀。由上可知：一幅图像具有较广的灰度级并且分布均匀，那么该图像具有较高的对比度和较大的灰度变化范围，这种图像可以提供更多的图像细节，更有利于进一步的图像处理。

直方图均衡化理论

图像直方图均衡化：均表示平均，即在图像中平均到每个灰度级，衡表示平衡、一样多，均衡表明直方图中每个灰度级上的直方图值一样大，直方图一样高，上图很形象地描述了均衡化，当然实际均衡化的效果不可能这样。
假设原始图像的灰度为$r$,目标灰度为$s$,因为灰度值是一个实数，则从数学的角度来看，实数的变化可以通过一个映射函数来实现，即:$$s=T(r),0\le r\le L-1$$该灰度变换函数应当满足以下条件：

• $0\le s=T(r)\le L-1$. 因为灰度变换后还是灰度，所以必须满足灰度值范围。
• $T(r)$在$0\le r\le L-1$上是一个严格单调递增函数。①在函数映射变化之后至少要求像素点的亮暗关系不变，比如说，$r_1$=30, $r_2$=20($r_1>r_2$)在灰度图中灰度$r_1$对应的像素点比灰度$r_2$对应的点亮，那么变换之后还应当满足$T(r_1)$对应的像素点比$T(r_2)$对应的像素点亮，即$T(r_1)>T(r_2)$，因此函数$T(r)$应该是增函数；②为什么需要严格的单调？在后面计算中需要用到$r=T^{-1}(s)$,该公式表示的是函数$T(r)$的反函数，计算反函数的前提就是要求$T(r)$单调，不单调则不存在反函数，同时也是保证$s$到$r$的反映射是一对一的(不严格单调则反映射存在一对多的情况)因此$T(r)$必须是单调函数。

上图左边为概率密度函数$p_r(r)$，在图像均衡化后最理想的效果就是图像的每个灰度级上的值一样大，再加上要求$p_s(s)$在$[0,L-1]$上的积分为1，所以最直观的$p_s(s)=\frac{1}{L-1}$
因为已知$r$和$s$存在的映射关系$T(r)$和$p_r(r)$，根据概率学知识就可以使用$p_r(r)$和$T(r)$来并表示$p_s(s)$，即：$$p_s(s)=p_r(r)|\frac{dr}{ds}|\text{\hspace{1em}(1)}$$上式中$r$为$s$的函数，即$r=T^{-1}(s)$
将$p_s(s)=\frac{1}{L-1}$代入式(1)可得灰度变换函数$T(r)$：$$\frac{1}{L-1}=p_r(r)|\frac{dr}{ds}|$$$$\frac{1}{L-1}ds=p_r(r)dr$$$${\int }_0^s\frac{1}{L-1}ds={\int }_0^r{p}_r(r)dr$$$$\frac{1}{L-1}s={\int }_0^r{p}_r(r)dr$$$$s=(L-1){\int }_0^r{p}_r(r)dr$$$$T(r)=(L-1){\int }_0^r{p}_r(r)dr$$$$T(r)=(L-1){\int }_0^r{p}_r(w)dw\text{\hspace{1em}(2)}$$

验证构造函数是否符合条件：
①${\int }_0^r{p}_r(w)dw$表示的是关于随机变量$r$的累计分布函数，其中$p_r(w)$D大于0，因此该积分项随着随机变量$r$的增大而增大，因此符合灰度变换函数的单调递增条件；
②当$r$=$L-1$时，$T(r)$增加到最大，为$L-1$，即${\int }_0^{L-1}{p}_r(w)dw$表示灰度直方图中灰度级0~$L-1$各自出现概率之和，为1，符合灰度变换函数的灰度范围条件。
由上述分析可知，$p_s(s)=\frac{1}{L-1}$是完全可行的概率密度函数，构造出来的灰度变换函数符合条件。由$p_s(s)=\frac{1}{L-1}$可知，它是一个均匀PDF，与$p_r(r)$的具体形式无关。
因为将数字图像的的各个灰度级视为一个离散随机变量，因此使用各个灰度级出现的概率(归一化后直方图的值)与求和来替代概率密度函数与积分，则等式(2)改写为：$$s_k=T(r_k)=(L-1)\sum _{j=0}^k{p}_r(r_j)=\frac{(L-1)}{MN}\sum _{j=0}^k{n}_j,k=0,1,2,\cdots ,L-1$$
等式中${\sum }_{j=0}^k{p}_r(r_j)$表示累计直方图值
图像直方图均衡化步骤：
①计算灰度值$r_k$出现的概率；
②计算累计直方图值，即${\sum }_{j=0}^k{p}_r(r_j)$；
③将累计直方图值代入公式计算；
④将计算结果取整；
⑤将新的灰度值放入与原灰度值对应的位置，即可生成均衡化后的图像。
实际例子计算：

如上表格所示，该图像为64*64的3比特数字图像，第一列为灰度级$r_k$，第二列为对应的灰度值个数$n_k$，第三列为直方图值$p_r(r_k)$，第四列为累计直方图值(归一化后的直方图和${\sum }_j^{L-1}{p}_r(r_k)$)，第五列为均衡化后的灰度值$s_k$。第六列为取整后的灰度值

直方图均衡化实现

直方图生成

• 生成灰度直方图

src = cv.imread('lenna.png', 0)#读取原始图像
src = cv.resize(src, (300, 300))#缩小图像尺寸为300*300

array = np.zeros([256])#新建空矩阵用于存放灰度值个数，初始化为0，256表示灰度级为256，
for i in range(src.shape[0]):
    for j in range(src.shape[1]):
        array[src[i, j]] += 1#遍历整张图像，若灰度值相同则就在矩阵中对应位置加1
# array = array /(300*300)#将数据矩阵中的每一个值都除以300*300，进行归一化处理

plt.title('histogram')
plt.plot(array, c='r')#用实线绘制
plt.hist(src.ravel(), 256)#利用函数绘制灰度图
plt.show()

• 实验结果:

图中红色部分(该部分使用实线绘制,故表现为连续,实际应是离散值)为自编代码实现,而蓝色部分为cv内部函数实现.

• 彩色直方图

src = cv.imread('lenna.png', 1)
array = np.zeros([256, 3])
for k in range(src.shape[2]):
    for i in range(src.shape[0]):
        for j in range(src.shape[1]):
                array[src[i, j, k],k] += 1

hist_b1 = cv.calcHist([src], [0], None, [256], [0, 255])#使用cv内部函数绘制B通道灰度图
hist_g1 = cv.calcHist([src], [1], None, [256], [0, 255])#使用cv内部函数绘制G通道灰度图
hist_r1 = cv.calcHist([src], [2], None, [256], [0, 255])#使用cv内部函数绘制R通道灰度图

plt.plot(hist_b1, c='b')
plt.plot(hist_g1, c='g')
plt.plot(hist_r1, c='r')

plt.plot(array[:, 0], 'y*')
plt.plot(array[:, 1], 'c*')
plt.plot(array[:, 2], 'k*')

plt.show()

效果对比:

(同样,实线部分并非真正直方图,只是为了更形象对比两种实现效果)

直方图均衡化

灰度直方图均衡化

def histogram_equalization_fnuc(src):
    src_array = np.zeros([256, 2])#建立零矩阵存储灰度值个数
    for i in range(src.shape[0]):
        for j in range(src.shape[1]):
            src_array[src[i, j], 0] = src[i, j]#矩阵的第一列用于存储灰度值
            src_array[src[i, j], 1] += 1#矩阵的第二列用于存储灰度值对应的数目
    src_array[:, 1] /= (300 * 300)#直方图归一化处理,则矩阵的第二列用于存储灰度值出现的概率

    accum = 0
    result_img = np.zeros_like(src)#建立一个矩阵用于存储均衡化后的数据

    for m in range(256):
        accum += src_array[m, 1]#计算累加直方图值(概率)
        temp = np.round(accum * 255)#将计算后的像素值取整(灰度图均衡化关键一步)
        for i in range(src.shape[0]):
            for j in range(src.shape[1]):
                if src[i, j] == src_array[m, 0]:#根据之前存储的灰度值,遍历原始图像,判断图像中与当前灰度相同的位置,
                                                # 若相同则在新的图像位置存储均衡化后的灰度值
                    result_img[i, j] = temp
    return result_img

主函数：

def main():
    src = cv.imread('lenna.png', 0)#读取灰度图像
    src = cv.resize(src, (300, 300))#缩小尺寸
    plt.hist(src.ravel(), 256, color='b')#查看原灰度图直方图

    result = histogram_equalization_fnuc(src)#对原灰度图进行均衡化
    plt.hist(result.ravel(), 256, color='r')#查看原灰度图均衡化后的直方图
    plt.show()
    cv.imshow("not_func", np.hstack([src, result]))#对比均衡化前后图像效果
    cv.waitKey(0)

实现效果：

上图中蓝色部分为原灰度图的直方图，红色部分为均衡化后的直方图，可以明显看出，均衡化之后图像的灰度级增多，几乎覆盖了8比特图全部灰度级，提高了图像对比度。

彩色直方图均衡化

def main():
    b, g, r = cv.split(src)#彩色图像均衡化需要分通道进行
    b = histogram_equalization_fnuc(b)
    g = histogram_equalization_fnuc(g)
    r = histogram_equalization_fnuc(r)
    result = cv.merge([b, g, r])#最后将均衡化后的三通道进行合并
    cv.imshow("Histogram Equalization", np.hstack([src, result]))#均衡化前后对比，
    cv.waitKey(0)
    cv.destroyAllWindows()

if __name__ == '__main__':
    main()

实现效果：

因为每个通道的灰度图都发生了均衡化，所以三个通道组合的彩色图颜色也发生了变化。

参考文献：[1]数字图像处理：第3版/(美)冈萨雷斯(Gonzalez,R.C),(美)伍兹(Woods, R.E.)著.阮秋琦等译
[2]直方图匹配的数学原理

[Until：2021年1月15日17时，转眼又是半个月……]