直方图均衡 Histogram Equalization

亮度直方图

在说明直方图均衡之前，先说说亮度直方图的概念。为了评估一个图像的色调转换，首先需要建立亮度直方图。亮度直方图就是图像中亮度分布的图表。在横轴上表示亮度值从黑色到白色；在竖轴上表示某一亮度所累积的像素数量。这里的亮度值指的是灰度等级，范围一般从 0 到 255，0 表示黑色，255 表示白色。

上面图片显示的是对比度差的两个例子。第一个图片，其亮度范围没有全部使用。在图表中可以看出，0 和 255 的位置上没有对应的亮度值，说明在图片中没有白色和黑色。第二个图片，亮度范围被全部使用，但是亮度聚集在某些峰值附近。所以，该图片大多数像素具有相同的亮度。

直方图均衡

有时候，因为摄像头传感器或者外界光线等原因，我们会得到像素值限定在某一范围内的图片，如上面的例子所示。但是，这样的图片并不是好的图片。好的图片应该是包含所有范围的亮度值。这时候就需要用到直方图均衡(Histogram Equlization)来处理这种情况，简单过程如下图所示：

来自维基百科

简单来说，直方图均衡化是使用图像直方图对对比度进行调整的图像处理方法。目的在于提高图像的全局对比度，使亮的地方更亮，暗的地方更暗。常被用于背景和前景都太亮或者太暗的图像，尤其是 X 光中骨骼的显示以及曝光过度或者曝光不足的图片的调整。

X 光中骨骼的显示，右边是均衡后的结果

实现过程

对于一个灰度图像 {x}，ni 表示为灰度级别为 i 的出现的次数。在图像中出现级别 i 的像素的概率为：

p_{x}(i) = p(x=i) = \frac{n_{i}}{n}, 0\leq i <L

L 是图像中灰度级别的总数（通常为256），n 是图像中的像素总数，px(i) 实际上是像素值 i 的图像直方图，归一化为 [0,1]。

直方图均衡化的处理依赖于累积概率函数（cdf）的使用。 cdf 是位于其域中的所有概率的累积和，数字图像的 cdf 定义如下：

cdf_{x}(i)=\sum_{j=0}^{i}p_{x}(j)

模拟图像的 cdf 定义如下：

cdf_{r}(r) = \int_{0}^{r}p_{r}(w)dw

该处理方法的思想是使用原图像的累积分布函数来转换像素值。下面使用模拟图像的 cdf 来证明为什么可以使用 cdf 来当转换函数。

首先，记住我们的目标是希望像素值的分布从左图变为右图的均匀分布。这样才能使得像素值在所有亮度范围中均匀分布，而不是集中在某一部分。

假设我们的转换函数为 cdf：

s = T(r) = cdf_{r}(r) = \int_{0}^{r}p_{r}(w)dw

通常，在概率统计中可以得到这样的定律：

P_{s}(s) = P_{r}(r)| \frac{dr}{ds}|

其中 Ps(s) 假设为新的概率分布，Pr(r) 为原来的概率分布。dr/ds 的关系可以从下面推导式中得出：

\frac{ds}{dr} = \frac{dT(r)}{dr}=(L-1)\frac{\int_{0}^{r}p_{r}(w)dw}{dr} = (L-1)·P_{r}(r)

将 dr/ds 的关系代入到下面式子中，可以发现 Ps(s) 的分布为均匀分布，符合我们的目标：

P_{s}(s) = P_{r}(r)| \frac{dr}{ds}| = \frac{1}{L-1}

因此，我们可以用下面的转换公式也就是 cdf，来得到新的像素值：

T(r) = (L-1) · cdf

通常，还可以使用下面这个转换函数：

T(r)= round(\frac{(cdf(r) - cdf_{min} )· (L-1) }{cdf_{max} - cdf_{min}})

上面就是灰度图片的直方图均衡实现过程。简单概括就是对像素进行点操作（point operator），将原始像素值映射为另一个范围。

分别将上面过程应用于 RGB 图像的红，绿，蓝通道，就可以对彩色图片进行直方图均衡处理。

但实际上，对彩色分量 RGB 分别做均衡化，会产生奇异的点，破坏图像的色彩平衡。一般采用的是用 HSL 和 HSV 色彩空间进行亮度的均衡即可。

应用 NumPy 实现

下面对一张整体较暗的图片进行处理。

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


img = cv2.imread('bay.jpg', 0)
hist, bins = np.histogram(img.flatten(), 256, [0, 256])

cdf = hist.cumsum()
cdf_normalized = cdf * hist.max() / cdf.max()

plt.plot(cdf_normalized, color='b')
plt.hist(img.flatten(), 256, [0, 256], color = 'r')
plt.xlim([0, 256])
plt.legend(('cdf', 'histogram'), loc='upper left')
plt.show()

原图片的直方图和 CDF 如下所示：

下面进行直方图均衡：

# calculate cdf
cdf_m = np.ma.masked_equal(cdf, 0)
cdf_m = (cdf_m - cdf_m.min())*255 / (cdf_m.max() - cdf_m.min())
cdf = np.ma.filled(cdf_m, 0).astype('uint8')

# mapping
img2 = cdf[img]
res = np.hstack((img, img2))
plt.figure(figsize=(10,10))
plt.imshow(res, cmap="gray")

结果如下所示，右边的是均衡后的结果：

下面计算均衡后的 CDF:

hist, bins = np.histogram(img2.flatten(), 256, [0, 256])

cdf = hist.cumsum()
cdf_normalized = cdf * hist.max() / cdf.max()

plt.plot(cdf_normalized, color='b')
plt.hist(img2.flatten(), 256, [0, 256], color = 'r')
plt.xlim([0, 256])
plt.ylim([0, 30000])
plt.legend(('cdf', 'histogram'), loc='upper left')
plt.show()

均衡后的直方图和 CDF 为：

可以看出均衡后的直方图比原先的亮度范围要宽，覆盖了整个亮度范围。

OpenCV 的实现过程

img = cv2.imread('bay.jpg', 0)
equ = cv2.equalizeHist(img)
res = np.hstack((img,equ)) 
cv2.imwrite('res.png', res)

plt.figure(figsize=(10,10))
plt.imshow(res, cmap="gray")

Skimage 的实现过程

from skimage import data, img_as_float
from skimage import exposure

img = plt.imread('bay.jpg')
img = img_as_float(img)
img_eq = exposure.equalize_hist(img)
res = np.hstack((img, img_eq))
cv2.imwrite('res.png', res)

plt.figure(figsize=(10,10))
plt.imshow(res, cmap="gray")

全部代码地址：histogram-equalization

参考

[1]. OpenCV - Histogram Equalization [2]. wiki - Histogram equalization [3]. Coursera - Image and Video Processing