直方图均衡化(II)

各位大家好，我是灿视，今天是直方图均衡的第二篇～

上一篇文章，我们主要是给大家看了下直方图均衡干了什么事情，并且直接给出了，针对离散型数据的直方图均衡化的公式。

今天，我们来看下，这里面的推导过程，为什么它的公式要这么做，以及改进的两类非全局直方图均衡化～。

1. 直方图均衡公式推导

在上一篇文章，我们了解到均衡化的目的是将原始图像的直方图变为均衡分布的的形式，将一非均匀灰度概率密度分布图像，通过寻求某种灰度变换，变成一幅具有均匀概率密度分布的目的图像。我们开始以直方图均衡化的推导公式入手：

假定：rrr代表灰度级，P(r)P(r)P(r)为概率密度函数。rrr值已经过归一化处理，灰度值范围在[0,1][0,1][0,1]之间。rrr与P(r)P(r)P(r)之间的关系如下左图所示。均衡化的目的是将上面的非均匀分布变成如下右图所示的均匀分布：

我们需要使用一种变换S=T(r)S=T(r)S=T(r)使直方图变平直，为使变换后的灰度仍保持从黑到白的单一变化顺序，且变换范围与原先一致，以避免整体变亮或变暗，此函数需要满足以下条件：

在0<=r<=10 <= r <= 10<=r<=1中，T(r)T(r)T(r)是单调递增函数，且0<=T(r)<=10 <= T(r) <= 10<=T(r)<=1；
反变换r=T−1(s),T−1(s)\mathrm{r}=T^{-1}(\mathrm{~s}), T^{-1}(\mathrm{~s})r=T−1( s),T−1( s)也为单调递增函数，且0<=s<=10 <= s <= 10<=s<=1。

具体的变换过程如下所示：

因为灰度变换不影响像素的位置分布，而且也不会增减像素数目，所以有如下的推导公式：
∫0rp(r)dr=∫0sp(s)ds=∫0s1⋅ds=s=T(r)T(r)=∫0rp(r)dr\begin{array}{l} \int_{0}^{r} p(r) d r=\int_{0}^{s} p(s) d s=\int_{0}^{s} 1 \cdot d s=s=T(r) \\ T(r)=\int_{0}^{r} p(r) d r \end{array} ∫0rp(r)dr=∫0sp(s)ds=∫0s1⋅ds=s=T(r)T(r)=∫0rp(r)dr

2. 离散灰度级

因此，对于离散型的图像：

设一幅图像的像素总数为nnn，分为LLL个灰度级，其中：

nkn_{k}nk：表示第K个灰度级出现的个数。

p(rk)=nk/np(r_{k})=n_{k}/np(rk)=nk/n：第KKK个灰度级出现的概率。

(0<=rk<=1，k=0,1,2,...,L−10<=r_{k}<=1，k=0,1,2,...,L-10<=rk<=1，k=0,1,2,...,L−1)公式如下：
Sk=T(rk)=∑j=0kPr(rj)=∑j=0knjnS_{k}=T\left(r_{k}\right)=\sum_{j=0}^{k} P_{r}\left(r_{j}\right)=\sum_{j=0}^{k} \frac{n_{j}}{n} Sk=T(rk)=j=0∑kPr(rj)=j=0∑knnj
因此，整个计算步骤如下：

求出图像中所包含的灰度级rkr_{k}rk都经过归一化处理，范围在[0,1][0,1][0,1]之间，也可以定在[0,L−1][0,L-1][0,L−1]之间。
统计各灰度级的像素数目nkn_{k}nk。
计算图像直方图。
计算变换函数，即Sk=T(rk)=∑j=0kPr(rj)=∑j=0knjnS_{k}=T\left(r_{k}\right)=\sum_{j=0}^{k} P_{r}\left(r_{j}\right)=\sum_{j=0}^{k} \frac{n_{j}}{n}Sk=T(rk)=∑j=0kPr(rj)=∑j=0knnj。
用变换函数计算映射后输出的灰度级sks_{k}sk。
统计映射后新的灰度级sks_{k}sk的像素数目nkn_{k}nk。
计算输出图像的直方图。

3. 举个例子

假设有64∗6464 * 6464∗64大小的图像，有888个灰度级，灰度分布如下所示：
rknkP(rk)r0=07900.19r1=1/710230.25r2=2/78500.21r3=3/76560.16r4=4/73290.08r5=5/72450.06r6=6/71220.03r7=1810.02\begin{array}{ccc} r_{k} & n_{k} & P\left(r_{k}\right) \\ r_{0}=0 & 790 & 0.19 \\ r_{1}=1 / 7 & 1023 & 0.25 \\ r_{2}=2 / 7 & 850 & 0.21 \\ r_{3}=3 / 7 & 656 & 0.16 \\ r_{4}=4 / 7 & 329 & 0.08 \\ r_{5}=5 / 7 & 245 & 0.06 \\ r_{6}=6 / 7 & 122 & 0.03 \\ r_{7}=1 & 81 & 0.02 \end{array} rkr0=0r1=1/7r2=2/7r3=3/7r4=4/7r5=5/7r6=6/7r7=1nk790102385065632924512281P(rk)0.190.250.210.160.080.060.030.02
由上，我们知道该图像的rkr_{k}rk，nkn_{k}nk和p(rk)p(r_{k})p(rk)。

其中：s0=T(r0)=0.19s_{0}=T(r_{0})=0.19s0=T(r0)=0.19

s1=T(r1)=0.19+0.25=0.44s_{1}=T(r_{1})=0.19 + 0.25=0.44s1=T(r1)=0.19+0.25=0.44

s2=T(r2)=0.19+0.25+0.21=0.65s_{2}=T(r_{2})=0.19 + 0.25 + 0.21=0.65s2=T(r2)=0.19+0.25+0.21=0.65

s3=T(r3)=0.19+0.25+0.21+0.16=0.81s_{3}=T(r_{3})=0.19 + 0.25 + 0.21 + 0.16=0.81s3=T(r3)=0.19+0.25+0.21+0.16=0.81

s4=T(r4)=0.19+0.25+0.21+0.16+0.08=0.89s_{4}=T(r_{4})=0.19 + 0.25 + 0.21 + 0.16 + 0.08=0.89s4=T(r4)=0.19+0.25+0.21+0.16+0.08=0.89

s5=T(r2)=0.19+0.25+0.21+0.16+0.08+0.06=0.95s_{5}=T(r_{2})=0.19 + 0.25 + 0.21 + 0.16 + 0.08 +0.06=0.95s5=T(r2)=0.19+0.25+0.21+0.16+0.08+0.06=0.95

s6=T(r2)=0.19+0.25+0.21+0.16+0.08+0.06+0.03=0.98s_{6}=T(r_{2})=0.19 + 0.25 + 0.21 + 0.16 + 0.08 +0.06+0.03=0.98s6=T(r2)=0.19+0.25+0.21+0.16+0.08+0.06+0.03=0.98

s7=T(r2)=0.19+0.25+0.21+0.16+0.08+0.06+0.03+0.02=1s_{7}=T(r_{2})=0.19 + 0.25 + 0.21 + 0.16 + 0.08 +0.06 +0.03 + 0.02=1s7=T(r2)=0.19+0.25+0.21+0.16+0.08+0.06+0.03+0.02=1

由于原图像的灰度级只有8级，变换之后的sks_{k}sk只能选择最接近的一个灰度级，因此需要对sks_{k}sk进行舍入处理，即上述步骤中的第五步，对每个sks_{k}sk将以1/71/71/7为量化单位进行舍入运算，结果如下：

s0=1/7,s1=3/7,s2=5/7,s3=6/7,s4=6/7,s5=1,s6=1,s7=1s_{0}=1 / 7, s_{1}=3 / 7, s_{2}=5 / 7, s_{3}=6 / 7, s_{4}=6 / 7, s_{5}=1, s_{6}=1, s_{7}=1s0=1/7,s1=3/7,s2=5/7,s3=6/7,s4=6/7,s5=1,s6=1,s7=1

根据舍入后的结果，我们可以得到均衡化后的灰度级仅有5个级别，分别是：

s0=1/7,s1=3/7,s2=5/7,s3=6/7,s4=1s_{0}=1 / 7, s_{1}=3 / 7, s_{2}=5 / 7, s_{3}=6 / 7, s_{4}=1s0=1/7,s1=3/7,s2=5/7,s3=6/7,s4=1

统计映射后新的灰度级 sks_{k}sk 的像素数目 nkn_{k}nk, 然后就可得到均衡化后的概率密度函数 p(sk)p\left(s_{k}\right)p(sk) 。
rknkP(rk)Sk计算Sk舍入SknskP(sk)r0=07900.190.191/7s07900.19r1=1/710230.250.443/7s110230.25r2=2/78500.210.655/7s28500.21r3=3/76560.160.816/7r4=4/73290.080.896/7s39850.24r5=5/72450.060.951r6=6/71220.030.981r7=1810.021.001s44480.11\begin{array}{cccccccc} r_{k} & n_{k} & P\left(r_{k}\right) & S_{k 计算} & S_{k舍入} & S_{k} & n_{s k} & P\left(s_{k}\right) \\ r_{0}=0 & 790 & 0.19 & 0.19 & 1 / 7 & s_{0} & 790 & 0.19 \\ r_{1}=1 / 7 & 1023 & 0.25 & 0.44 & 3 / 7 & s_{1} & 1023 & 0.25 \\ r_{2}=2 / 7 & 850 & 0.21 & 0.65 & 5 / 7 & s_{2} & 850 & 0.21 \\ r_{3}=3 / 7 & 656 & 0.16 & 0.81 & 6 / 7 & & & \\ r_{4}=4 / 7 & 329 & 0.08 & 0.89 & 6 / 7 & s_{3} & 985 & 0.24 \\ r_{5}=5 / 7 & 245 & 0.06 & 0.95 & 1 & & & \\ r_{6}=6 / 7 & 122 & 0.03 & 0.98 & 1 & & & \\ r_{7}=1 & 81 & 0.02 & 1.00 & 1 & s_{4}& 448 & 0.11 \end{array} rkr0=0r1=1/7r2=2/7r3=3/7r4=4/7r5=5/7r6=6/7r7=1nk790102385065632924512281P(rk)0.190.250.210.160.080.060.030.02Sk计算0.190.440.650.810.890.950.981.00Sk舍入1/73/75/76/76/7111Sks0s1s2s3s4nsk7901023850985448P(sk)0.190.250.210.240.11
因此，直方图均衡化实质上是减少图像的灰度级来加大对比度，图像经均衡化处理之后，图像变得清晰，直方图中每个像素点的灰度级减少，但分布更加均匀，对比度更高。

如上文所示的直方图均衡之后的效果：

如上图，直方图均衡化操作是对全局进行均衡化，直方图覆盖的范围太大，图片变得太亮了，反而会丢失的一些图像信息。

4. 自适应直方图均衡

在前面介绍的直方图均衡化中，是直接对全局图像进行均衡化，是GlobalGlobalGlobal HistogramHistogramHistogram EqualizationEqualizationEqualization，而没有考虑到局部图像区域(LocalLocalLocal RegionRegionRegion)，自适应过程就是在均衡化的过程中只利用局部区域窗口内的直方图分布来构建映射函数fff 。

最早的自适应直方图均衡(AHEAHEAHE(AdaptiveAdaptiveAdaptive HisogramHisogramHisogram EqualizationEqualizationEqualization))方案如下所示：

对A图像每个像素点进行遍历，用像素点周围W∗WW * WW∗W的窗口进行计算直方图变换的分布函数，然后对该像素点进行映射。但首先需要对窗口大小为W∗WW * WW∗W的图像计算分布函数，复杂度为O(W∗W+L)O(W*W + L)O(W∗W+L) ，因此，对整张图像遍历一遍下来，复杂度为O(M∗N∗(W∗W+L))O(M*N*(W*W+L))O(M∗N∗(W∗W+L))，复杂度太高了。

针对上面高复杂度的问题，给出了一系列改进的方案。在之前使用W∗WW * WW∗W的窗口来计算直方图分布函数，之后不仅对图像的一个像素点进行变换，而是将一些列(如，4， 9， 16…)像素点进行变化。若是利用64∗6464*6464∗64的窗口来进行计算，此时采用32∗3232 * 3232∗32的中心区域像素点进行变换，这样的话，速度就是Naive AHE的32∗3232*3232∗32倍。相当于在做图像分快的同时，做了直方图均衡化。

但是，这样的操作，会产生一种马赛克现象，这是因为当W∗WW * WW∗W窗口内像素点近似一样时，即直方图只有一个灰度级，那么这样得到的分布函数就会是一种“阶跃”的曲线，使变换后图像过度增强，一些噪声就可能被过度放大。

5. 限制对比度自适应直方图均衡化(CLAHE)

我们可以先将图像分块，每块计算一个直方图分布函数，即如下的黑丝实线边框的小块。之后，分别计算四个窗口直方图的分布函数对蓝色像素点的映射值，ful(D),fur(D),fdl(D),fdr(D)f_{u l}(D), f_{u r}(D), f_{d l}(D), f_{d r}(D)ful(D),fur(D),fdl(D),fdr(D)。再对其进行双线性插值，其公式如下所示：

f(D)=(1−Δy)((1−Δx)ful(D)+Δxfbl(D))+Δy((1−Δx)fur(D)+Δxfbr(D))f(D)=(1-\Delta y)\left((1-\Delta x) f_{u l}(D)+\Delta x f_{b l}(D)\right)+\Delta y\left((1-\Delta x) f_{u r}(D)+\Delta x f_{b r}(D)\right)f(D)=(1−Δy)((1−Δx)ful(D)+Δxfbl(D))+Δy((1−Δx)fur(D)+Δxfbr(D))
Δx,Δy\Delta x, \Delta yΔx,Δy 是蓝色像素点相对于左上角窗口中心，即左上角黑色像素点的距离与窗口大小的比值。

对于图中红色像素点，只使用其最近的窗口的CDF进行映射；对于图中绿色像素点，采用邻近的两个窗口的CDF映射值进行线性插值。

CLAHECLAHECLAHE涉及到窗口大小、影响区域大小、直方图阈值三个部分的参数，可以通过调参更好地去控制对比度增强的结果。

图像被划分后的小矩形框被称为 “tiles” (在 OpenCV 中，tileSize 默认为8x8)，在这每个小矩形框内，直方图也会被限制在一个很小的范围内（除非图像中有噪声）。

如果小区域中存在噪声，那么该噪声就会被放大，为了避免这一点，使用了对比度限制 (在 OpenCV 中，clipLimit 默认为 2.0)。

img = cv2.imread('boy.png', 0)# create a CLAHE object
clahe = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(8,8))
cl1 = clahe.apply(img)
res = np.hstack((img, cl1))
cv2.imshow("dst", res)
cv2.waitKey(0

6.自适应局部区域伸展(Local Region Stretch)直方图均衡化

上面提到的AHEAHEAHE和CLAHECLAHECLAHE都是基于块状区域进行直方图均衡化的，但是能不能根据灰度级区域近似的区域进行均衡化呢？比如对图像中灰度级[min,max][min, max][min,max]范围里面的所有像素点进行均衡化，使得像素点的直方图尽量在[min,max][min, max][min,max]上均匀分布。

因此，我们需要统计图像直方图，按照灰度级划分为三个灰度区间，使得三个区间的像素点数量近似相等，这样就分别在[0,level1)[0, level1)[0,level1),$ [level1, level2),,, [level2, 255]$三个灰度区间做直方图均衡化，最后合并。

如下图所示：

7. 小结

这篇主要是接着上一篇文章，推导了下直方图均衡化的公式。再引出自适应直方图均衡化(AHEAHEAHE) 以及限制对比度自适应直方图均衡化(CLAHECLAHECLAHE) 等直方图均衡化算法。

8. 参考

https://zhuanlan.zhihu.com/p/44918476
https://blog.csdn.net/qq_42593220/article/details/103227609
https://zhuanlan.zhihu.com/p/98541241