直方图均衡

直方图均衡与精确直方图均衡

直方图均衡化是一种图像处理技术，其基本原理是通过重分配图像中每个像素的灰度级别来扩展其灰度差异，使图像的灰度级别分布更均匀，整个图像具有更广泛的灰度范围，从而提高图像的对比度。

这里就简单介绍一下常用的灰度直方图均衡方法，并拓展至彩色直方图均衡、精确直方图均衡及其实现步骤~

灰度直方图均衡

在进行灰度直方图均衡之前，需要先计算图像的灰度直方图。

灰度直方图是一个显示图像中每个灰度值出现次数的图表，对于灰度图像来说，每个像素点是由0-255范围的灰度值表示的。计算每个灰度值在图像中出现的次数，并将结果绘制为条形图，就形成了灰度直方图。就好像下面这张图：

他的灰度直方图是这样的：

这张图是不是看起来很眼熟？如果经常摄影或者用PS的话，应该会看过不少次这种类型的图片。点开手机相册中随便一张图，点一下右上角的i，就可以看到图片的详细信息，弹出的信息里面就有各个颜色通道的直方图。回到这张图上来，原图的画面本身较暗，对比度不高，所以其像素的灰度值均集中在100以下，且排列非常紧密。

那么怎么实现直方图均衡呢？其实很简单，只需要让直方图上每个灰度对应的一列像素往右边挪一挪就可以了，这也是直方图均衡最基本的思想。最简单粗暴的实现方法就是把每个灰度值依次映射到另一个高一点的灰度值上，将对应像素的灰度值改成新的灰度值，就大功告成了。

实现步骤

首先读取图像，并将彩色图像转换为灰度图像。随后使用cv2.calcHist计算灰度图像的直方图。由于这里只需要计算灰度图像的直方图，因此通道索引为[0]，不使用掩码，直方图的大小为256，范围在0到255之间。

import numpy as np
import cv2
import matplotlib.pyplot as plt

# 选择彩色图像并加载，转化为灰度图像
image_path = 'image.jpg'
image = cv2.imread(image_path)
image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)

# 计算灰度直方图
hist = cv2.calcHist([image], [0], None, [256], [0, 256])

之后使用cumsum函数对直方图进行累积求和，得到累积分布函数（CDF），并对其进行归一化处理，将CDF的值乘以255，再除以CDF的最后一个值，使CDF映射到0-255的范围内。

cdf = hist.cumsum()
cdf_normalized = 255 * cdf / cdf[-1]  # 归一化处理

然后使用np.interp函数实现了线性插值，根据对应的灰度级别在CDF中的映射值进行插值，将原始图像中每个像素的灰度值映射为新的灰度值。

# 使用累计分布函数进行均衡化
transformed_image_array = np.interp(image.ravel(), range(256), cdf_normalized).reshape(image.shape)

最后使用plt.hist函数绘制原始灰度直方图和均衡后的灰度直方图，并展示处理前后的两幅图片。

# 绘制原始灰度直方图和均衡后的灰度直方图
plt.subplot(1,2,1)
plt.hist(image.flatten(), 256, [0, 256], color='gray')
plt.title('Original Histogram')

plt.subplot(1,2,2)
plt.hist(transformed_image_array.flatten(), 256, [0, 256], color='gray')
plt.title('Transformed Histogram')

plt.tight_layout()
plt.show()

# 显示原始图像和处理后的图像
plt.subplot(1,2,1)
plt.axis('off')
plt.imshow(image, cmap='gray')
plt.title('Original Image')

plt.subplot(1,2,2)
plt.imshow(transformed_image_array, cmap='gray')
plt.axis('off')
plt.title('Transformed Image')

plt.tight_layout()
plt.show()

（实际上真正用到的只有中间四句话...甚至还没plt的四分之一多...）

实现效果

图片的暗部细节明显清晰了不少，对比度明显提升。直方图看起来就更清楚了，一下就能看懂均衡的方法（非常的简单粗暴啊家人们，效果非常的直观哈）

彩色直方图均衡

灰度直方图会了，那拓展到彩色直方图不也轻轻松松？无非就是把原先的一个单通道，拓展成了RGB三个通道，也就是说只需要对三个通道分别进行之前的计算，不就实现了。

实现步骤

（进行一个很快的实现）这里只多了两个步骤，一个是分离通道，也就是

# 对每个颜色通道进行直方图均衡
red_channel = image[:, :, 0]
green_channel = image[:, :, 1]
blue_channel = image[:, :, 2]

另一个就是合并通道，也就是

# 合并处理后的颜色通道
transformed_image_array = cv2.merge([red_transformed, green_transformed, blue_transformed]).astype('int')

red_transformed, green_transformed, blue_transformed 这三个值就是经过上文的方法处理过的RGB通道的像素，这里我就不再重复了（就问你快不快吧）

实现效果

目前仅常见的彩色图像直方图均衡方法就有四种:基于颜色分量的独立均衡法、基于三维联合概率的均衡法、基于HSI颜色空间的I分量或者SI均衡法和基于三维彩色直方图均衡化的彩色图像增强算法。这里说到的方法就是基于颜色分量的独立均衡法，感觉应该是最简单明了的方法。这里的方法仅用于示范，实际上更好更常用的应该算是基于HSV颜色空间的均衡，直接使用RGB的话可能会有失真等问题存在。

虽然前面一直在说“简简单单”“轻轻松松”什么的，实际上，涉及到彩色直方图时还是有不少问题的。不同于灰度直方图，通道数量越多，图片的信息量也就越大。如何能使得图片的细节更还原，视觉表现更好，或者更符合特定任务处理的要求，还是存在不少需要攻克的难题的。从我个人浅薄的见解来看，这些方法啊、效果啊，还是随着科技、人文的发展需要来不断改进的。

精确直方图均衡

不知道各位在做的时候有没有觉得，之前说的这些均衡方式很奇怪？如果只是进行插值、线性的映射，就好像只是把直方图从紧密的离散值，拉成了稀疏一点的离散值，本质上好像还没什么变化。那有没有一种办法，可以在不大幅改变图片的视觉效果的情况下，让直方图呈现出更加均匀的分布状态，或者换句话说，让直方图变化成指定的形状？

答案当然是肯定的。不过如果直接利用上文提到的方法，往往并不能产生精确形状的直方图。因为在进行均衡化的时候，有一些像素的数量发生了变化，或者在使用其他方法时，可能会因为像素值的四舍五入，使得一些像素被丢弃。这两种情况导致直方图在y轴发生了一定的偏移。这时候就需要精确直方图均衡，或者说叫精确直方图匹配。

这里首先推荐一篇IEEE上的论文Exact histogram specification ，就是专门说这个做法的。虽然已经是快二十年前的论文了，但是讲的还是蛮好的，使用的方法也很有意思，这里我就只聊聊他所提到的方法。

原理

首先对于输入的N×M大小的图像f，假设图像中每个像素的灰度级别为L，以及待实现的直方图H=h0,h1,...,hL−1（其中hi表示灰度级别 i 的像素数量）。

定义一个严格排序的像素序列，使得像素的排序关系与灰度级别的大小有关。由于图像的像素数量往往远大于灰度值，也就是说在排序中会出现很多相同灰度值的像素，这在后续像素值重新映射图像时可能会导致像素被赋予了不合适的值，使其与邻域像素不匹配，最后的图像结果就会出现问题。

为了解决这种问题，论文中提到了一个非常重要的做法，对图像进行卷积操作，即将某点像素和周围相邻的像素与对应的算子求卷积和，并将卷积和结果作为中心像素的灰度值进行记录。文中列出来6个算子，分别为

$$ ϕ1=[1]ϕ2=15[010111010]ϕ3=19[111111111]ϕ4=113[0010001110111110111000100]ϕ5=121[0111011111111111111101110]ϕ6=125[1111111111111111111111111]

$$ 经过卷积后，可以依次得到1幅原图像和5幅处理后的图像。以原图像为第一顺位，如果有相同的灰度值则根据其他图像来排列其先后顺序，从第二幅图像开始寻找对应位置的灰度值进行排序，如果仍然相同则继续在下一幅图像中寻找。

通过这些操作，我们能得到相对合适的排列顺序，并根据对应的索引值，按照需要匹配的直方图将像素灰度值进行映射。将排序后的像素序列分割成L组，每组包含hj个像素，其中j为组的索引。对于每一组像素，将其灰度级别设置为j，这样就可以得到最终的图像。

实现步骤

首先命名好对应的算子（由于第一个算子的计算结果就是图像本身，这里就不再重复计算）

matrix1 = np.array([[0,1,0],
                    [1,1,1],
                    [0,1,0]])/5
matrix2 = np.array([[1,1,1],
                    [1,1,1],
                    [1,1,1]])/9
matrix3 = np.array([[0,0,1,0,0],
                    [0,1,1,1,0],
                    [1,1,1,1,1],
                    [0,1,1,1,0],
                    [0,0,1,0,0]])/13
matrix4 = np.array([[0,1,1,1,0],
                    [1,1,1,1,1],
                    [1,1,1,1,1],
                    [1,1,1,1,1],
                    [0,1,1,1,0]])/21
matrix5 = np.array([[1,1,1,1,1],
                    [1,1,1,1,1],
                    [1,1,1,1,1],
                    [1,1,1,1,1],
                    [1,1,1,1,1]])/25

随后将图像中每一个像素及周围相应数量的像素取出为一个矩阵，与指定算子进行卷积，并将卷积和作为中心像素的灰度值，最后返回由该像素值所组成的矩阵。为防止卷积时超出图像范围，提前将图片四周扩展为[1/2算子宽度]个灰度值为0的像素。

def mul_Image(matrix):
    result = np.zeros_like(image)
    num = int(len(matrix)/2)
    expanded_size = (image.shape[0]+num*2, image.shape[1]+num*2)

    eximg = np.zeros(expanded_size, dtype=np.uint8)
    eximg[num:expanded_size[0]-num, num:expanded_size[1]-num] = image

    for i in range(num, image.shape[0]-num):
        for j in range(num, image.shape[1]-num):
            mat = image[i-num:i+num+1, j-num:j+num+1]
            res = np.multiply(mat, matrix)
            result[i,j] = int(np.sum(res))
    
    return result

对所有算子进行上述操作后，就可以得到1幅原图像和5幅处理后的图像，为了方便后续计算，将图像展开成一维数组。在这里还要用zip将所有数组压缩成一个元组并转为列表形式。

a = image.reshape(1,-1)[0]
a1 = mul_Image(matrix1).reshape(1,-1)[0]
a2 = mul_Image(matrix2).reshape(1,-1)[0]
a3 = mul_Image(matrix3).reshape(1,-1)[0]
a4 = mul_Image(matrix4).reshape(1,-1)[0]
a5 = mul_Image(matrix5).reshape(1,-1)[0]
arr = list(zip(a, a1, a2, a3, a4, a5))

之后用enumerate() 函数将列表转化为包含元素和其所在索引值的元组列表，再用sorted()进行排序。这里使用lambda函数指定排序的关键字，将六幅图像的灰度值按照上面的方法依次比较并排序。排序后就可以从列表中提取出原图像展开的一维数组中每个像素对应索引值。

sorted_arr = sorted(enumerate(arr), key=lambda x: (x[1][0], x[1][1], x[1][2], x[1][3], x[1][4], x[1][5]))
sorted_indices = [idx for idx, _ in sorted_arr]

利用索引值将一个空的一位列表中的位置进行灰度的映射。这里我将输出的灰度直方图指定为直线（每个灰度的像素数量相等）。将索引值序列平均分割成256组，每组的灰度值依次为0至255，并将索引值对应的像素进行修改，并reshape成图像原来的大小，就得到了最终的输出图像。

h,w = image.shape
res = np.zeros(h*w)
num=0
maxlen = int(h*w/256)
pix = 0
for index in sorted_indices:
    num+=1
    if num>maxlen:
        num=0
        pix+=1
    res[index]=pix
    
res = res.reshape(h,w)

这里的res就是得到的图像结果。后面的图片展示这里就不再重复了，直接plt就可以了。

实现效果

处理前后的图片效果：

处理前后的图像的灰度直方图：

可以很直接地看出，这里的输出图像和第一段实现的效果相差不大，都只是感觉暗部提亮了一点，对比度高了一点，但灰度直方图已经匹配为目标直方图（常数函数）。如果想要实现对任意形状的直方图进行匹配，那么只需要将最后映射时的直方图纵坐标的值与对应位置上的函数值（或其他图像直方图对应灰度上的数量值）就可以了。