基于python的对比度增强(线性变换、直方图正规化、直方图均衡化、CLAHE)

线性变换

假设输入图像为I，宽为W，高为H，输出图像为O，图像的线性变换可以用以下公式定义:
O(r,c)=a×I(r,c)+b,0≤r<H,0≤c<WO(r, c) = a × I(r, c) + b, 0 ≤ r < H, 0 ≤ c < W O(r,c)=a×I(r,c)+b,0≤r<H,0≤c<W
当a=1，b=0时，O为I的一个副本；如果a>1，则输出图像O的对比度比I有所增大；如果0<a<1，则O的对比度比I有所减小。而b值的改变，影响的是输出图像的亮度，当b>0时，亮度增加；当b<0时，亮度减小

import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt# 统计灰度直方图并绘制
def calGrayHist(I):h, w = I.shapegrayHist = np.zeros(256, np.uint64)for i in range(h):for j in range(w):grayHist[I[i][j]] += 1plt.plot(grayHist)plt.xlabel("gray label")plt.ylabel("number of pixels")img_gray = cv.imread("./001_07.bmp", 0)
out = 2.0 * img_gray + 5
# 进行数据截断，大于255的像素值截断为255
out[out > 255] = 255
out = out.astype(np.uint8)plt.figure(figsize=(20,5))
plt.subplot(141)
plt.imshow(img_gray, cmap = "gray")
plt.subplot(142)
calGrayHist(img_gray)
plt.subplot(143)
plt.imshow(out, cmap = "gray")
plt.subplot(144)
calGrayHist(out)
# plt.savefig("./showimg.jpg")
plt.show()

线性变换的参数需要根据不同的应用及图像自身的信息进行合理的选择，可能需要进行多次测试，所以选择合适的参数是相当麻烦的。直方图正规化就是基于当前图像情况自动选取a和b的值的方法

直方图正规化

# 直方图正规化
# 计算原图中出现的最小灰度级和最大灰度级
Imin, Imax = cv.minMaxLoc(img_gray)[:2]
'''
# 也可以使用numpy计算
Imax = np.max(img_gray)
Imin = np.min(img_gray)
'''
Omin, Omax = 0, 255
a = float(Omax - Omin) / (Imax - Imin)
b = Omin - a * Imin
out = a * img_gray + b
out = out.astype(np.uint8)

代码中计算原图中出现的最小灰度级和最大灰度级可以使用OpenCV提供的函数

minVal, maxVal, minLoc, maxLoc = cv.minMaxLoc(src[, mask])

返回值分别为：最小值，最大值，最小值的位置索引，最大值的位置索引。

正规化函数normalize: dst=cv.normalize(src, dst[, alpha[, beta[, norm_type[, dtype[, mask]]]]])

使用函数normalize对图像进行对比度增强时，经常令参数norm_type=NORM_MINMAX，此函数内部实现和上边讲的计算方法是相同的，参数alpha相当于O_max，参数beta相当于O_min。注意，使用normalize可以处理多通道矩阵，分别对每一个通道进行正规化操作。使用该函数的代码如下，实现结果和上边是相同的

out = np.zeros(img_gray.shape, np.uint8)
cv.normalize(img_gray, out, 255, 0, cv.NORM_MINMAX, cv.CV_8U) # 调库
cv.imshow("img", img)
cv.imshow("out", out)
cv.waitKey(0)
cv.destroyAllWindows()

直方图均衡化(Histogram Equalization)

直方图均衡化是指: 利用图像直方图对对比度进行调整的方法
直方图均衡化通常用来增加许多图像的局部对比度，尤其是当图像的有用数据的对比度相当接近的时候
为讨论方便起见，以 r 和 s 分别表示归一化了的原图像灰度和经直方图均衡化后的图像灰度。当 r = s = 0 时，表示黑色；当 r = s = 1 时，表示白色；当 r, s ∈ (0, 1) 时，表示像素灰度在黑白之间变化
在 [0,1] 区间内的任何一个 r ，经变换函数 T(r)T(r)T(r) 都可以产生一个对应的 s
s=T(r)s = T(r)s=T(r)
式中，T(r)T(r)T(r) 应当满足以下两个条件:

在 0 ≤ r ≤ 1 内，T(r)T(r)T(r) 为单调递增函数；(此条件保证了均衡化后图像的灰度级从黑到白的次序不变)
在 0 ≤ r ≤ 1 内有 0≤T(r)≤10 ≤ T(r) ≤ 10≤T(r)≤1。(此条件保证了均衡化后图像的像素灰度值在允许的范围内)

如果已知随机变量 r 的概率密度是 pr(r)p_{r}(r)pr(r)，而随机变量 s 是 r 的函数，则 s 的概率密度 ps(s)p_{s}(s)ps(s) 可以由 pr(r)p_{r}(r)pr(r) 求出

则有: ps(s)ds=pr(r)drp_{s}(s)ds = p_{r}(r)drps(s)ds=pr(r)dr
又有：从人眼视觉特性来考虑，一幅图像的灰度直方图如果是均匀分布的，那么该图像看上去效果比较好 (参考冈萨雷斯数字图像处理3.3节)。因此要做直方图均衡化，这里的 ps(s)p_{s}(s)ps(s) 应当是均匀分布的概率密度函数

由概率论知识可知，对于区间 [a,b] 上的均匀分布，其概率密度函数等于 1b−a\frac{1}{b-a}b−a1。如果原图像没有进行归一化，即 r∈[0,L−1]r \in [0, L-1]r∈[0,L−1], 那么 ps(s)=1(L−1)−0=1L−1p_{s}(s) = \frac{1}{(L-1)-0} = \frac{1}{L-1}ps(s)=(L−1)−01=L−11，归一化之后 r∈[0,1]r \in [0, 1]r∈[0,1]，所以这里的 ps(s)=11−0=1p_{s}(s) = \frac{1}{1-0} = 1ps(s)=1−01=1

将ps(s)=1p_{s}(s) = 1ps(s)=1代入上式有: ds=pr(r)drds = p_{r}(r)drds=pr(r)dr，
两边积分得: s=T(r)=∫0rpr(r)drs = T(r) = \int_{0}^{r}p_{r}(r)drs=T(r)=∫0rpr(r)dr
由于数字图像灰度级是离散的，所以推广到离散情况，即用频率代替概率，有:
sk=T(rk)=∑i=0kpr(ri)=∑i=0kniNs_{k} = T(r_{k}) = \sum_{i=0}^{k}p_{r}(r_{i}) = \sum_{i=0}^{k}\frac{n_{i}}{N}sk=T(rk)=i=0∑kpr(ri)=i=0∑kNni
式中，0⩽rk⩽1,k=0,1,2,...,L−10 \leqslant r_{k} \leqslant 1, k = 0, 1, 2, ..., L-10⩽rk⩽1,k=0,1,2,...,L−1 (注: 这里的 rk=kL−1r_{k} = \frac{k}{L-1}rk=L−1k，表示归一化后的灰度级；k表示归一化前的灰度级)，而且需要注意的是，这里的 sks_{k}sk 也是归一化后的灰度级，其值在 0 到 1 之间；有时需要将其乘以L-1再取整，使其灰度级范围在 0 到 L-1 之间，与原图像一致

# 全局直方图均衡化
def equalHist(img, z_max = 255): # z_max = L-1 = 255# 灰度图像矩阵的高、宽H, W = img.shape# S is the total of pixelsS = H * Wout = np.zeros(img.shape)sum_h = 0for i in range(256):ind = np.where(img == i)sum_h += len(img[ind])z_prime = z_max / S * sum_hout[ind] = z_primeout = out.astype(np.uint8)return out

从结果图的直方图分布来看，并不是理想的一条水平线，原因是我们是在连续分布上推导得到的转换函数，但作用于离散分布，且不允许出现非整数形式的新结果
不过，整体看来该转换函数使得新的分布具有展开直方图分布的效果，增强了整体的对比度

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

关于CLAHE原理讲解的大多数文章讲的都很泛，就像蜻蜓点水般随便扯一下就直接调库实现了，所以真要弄懂的话还是要看源码
http://www.realtimerendering.com/resources/GraphicsGems/gemsiv/clahe.c

上述HE算法得到的结果存在一些问题: 1) 部分区域由于对比度增强过大，成为噪点；2) 一些区域调整后变得更暗/更亮，丢失细节信息

针对问题1)，有人提出了CLHE (即 Contrast Limited HE对比度限制的HE算法)
对统计得到的直方图进行裁剪，使其幅值低于某个上限，将裁剪掉的部分均匀分布在整个灰度区间，以确保直方图总面积不变

可以看到，这时直方图又会整体上升了一个高度，会超过我们设置的上限
其实在具体实现的时候有很多解决方法，可以多重复几次裁剪过程，使得上升的部分变得微不足道，或是用另一种常用的方法:

设裁剪值为ClipThreshold，求直方图中高于该阈值的部分的总和totalExcess，此时如果把totalExcess均分给所有灰度级，那么最终的直方图整体会上升一个高度，这个高度是averageIncrease = totalExcess / N，以upper = ClipThreshold - averageIncrease为界限对直方图进行如下处理:
1.若幅值高于ClipThreshold，直接设置为ClipThreshold
2.若幅值处于upper和ClipThreshold之间，从totalExcess中取出ClipThreshold - 幅值，将幅值填补成ClipThreshold
3.若幅值低于upper，直接加上averageIncrease即可
实现完上面的操作后，还有一些剩余的totalExcess没有被分出去，这个剩余是来自于1.和2.，这时需要再把这些均匀分给那些目前幅值依旧低于ClipThreshold的灰度值

'''This function performs clipping of the histogram and redistribution of bins.The histogram is clipped and the number of excess pixels is counted. Afterwardsthe excess pixels are equally redistributed across the whole histogram (providing the bin count is smaller than the cliplimit).
'''
def ClipHistogram(pHistogram, ClipThreshold):totalExcess = 0# 累计超出阈值的部分for i in range(256):if pHistogram[i] > ClipThreshold:totalExcess += (pHistogram[i] - ClipThreshold)averageIncrease = totalExcess // 256upper = ClipThreshold - averageIncrease# 修剪直方图并重新分配数值for i in range(256):if pHistogram[i] >= ClipThreshold: # 幅值高于ClipThresholdpHistogram[i] = ClipThresholdelse:if pHistogram[i] > upper: # 幅值介于upper和ClipThreshold之间totalExcess -= (ClipThreshold - pHistogram[i])pHistogram[i] = ClipThresholdelse: # 幅值小于uppertotalExcess -= averageIncreasepHistogram[i] += averageIncreasepStartBin = 0while totalExcess > 0:# 上述过程后totalExcess仍会有剩余未分配，设置步长将剩余的进行均分step_size = int(256 / totalExcess)if step_size < 1: # step_size至少为1step_size = 1for i in range(pStartBin, 256, step_size):if totalExcess == 0:breakif pHistogram[i] < ClipThreshold:pHistogram[i] += 1totalExcess -= 1 # 减少totalExcesspStartBin += 1 # 在其他位置重新开始分配return pHistogram

针对问题2)，有人提出AHE (Adaptive HE自适应直方图均衡化算法)，基本思想是: 对原图中每个像素，计算其周围一个邻域内的直方图，并使用HE映射得到新的像素值。为了能够处理图像边缘的像素，一般先对原始图像做镜像扩边处理
不过这样做计算量偏大，为了减少计算量，一般先把原始图像分块，分别对每个子区域进行HE变换。但是这样又产生了新的问题，子区域相接处像素值分布不连续，产生很强的分割线

为了解决这个问题，提出了优化方案双线性插值的AHE，然后在这个基础上，再使用CLHE的限制对比度的思想，就变成了我们熟知的CLAHE算法

CLAHE算法的步骤:

图像分块，以块为单位，统计子块的直方图
使用限制对比度方法进行修剪直方图，最后再直方图均衡化
像素点灰度值重构，使用的是双线性插值

【双线性插值】
对图像分块，然后对每一块求直方图均衡化，但不是把直方图均衡化的结果全部赋给这个子块，因为这样会造成块跟块之间像素值不连续
那像素点的灰度级该如何更新? 根据像素点所在位置分成三类
在粉色区域内的点的灰度就等于直方图均衡化的结果，即直接由映射函数计算得到
在绿色区域内的点的灰度由相邻两个子图映射函数插值得到
在其他区域内的所有点的灰度由相邻四个子图映射函数双线性插值得到

如下图所示便是双线性插值的情况，O点的灰度值为ror_oro，对其周围最近的四个子块的灰度映射函数进行双线性插值的结果为:
so=(1−x)(1−y)TA(ro)+x(1−y)TB(ro)+(1−x)yTC(ro)+xyTD(ro)s_o = (1 - x)(1 - y)T_A(r_o) + x(1 - y)T_B(r_o) + (1 - x)yT_C(r_o) + xyT_D(r_o)so=(1−x)(1−y)TA(ro)+x(1−y)TB(ro)+(1−x)yTC(ro)+xyTD(ro)

【线性插值代码实现中的细节】
举个例子:
把一张图片分成4×3个子图，每个格子为4×4大小，图中的P、Q两点位于边界处，所以需要使用线性插值得到这两个点的灰度级，线性插值第一步找到该点周围相邻的两个子图，那P点和Q点的相邻子图分别是什么？假设P点坐标为(11, 2)、Q点坐标为(9, 2)

可以看到，虽然P点和Q点在同一个子图中，但相邻子图是不同的，Q点的灰度级是由TA(r)T_A(r)TA(r)和TB(r)T_B(r)TB(r)这两个映射函数得到，也就是说Q点相邻的两个子图是第3个和第6个子图(从0开始编号)，而P点的灰度级是由TB(r)T_B(r)TB(r)和TC(r)T_C(r)TC(r)这两个映射函数得到，是第6个子图和第9个子图
代码中如何判定？

Q = (9, 2)
P = (11, 2)
tileGridSize = (4, 3) # 将图片分成4×3份
height_block, width_block = 4, 4 # 每个子图的大小Q_numi = int((Q[0] - height_block / 2) / height_block)
Q_num1 = Q_numi * tileGridSize[1]
Q_num2 = Q_num1 + tileGridSize[1]P_numi = int((P[0] - height_block / 2) / height_block)
P_num1 = P_numi * tileGridSize[1]
P_num2 = P_num1 + tileGridSize[1]print(Q_num1, Q_num2)
print(P_num1, P_num2)

运行结果为:

3 6
6 9

Q点的灰度级公式就为:
sQ=(1−yQ)TA(rQ)+yQTB(rQ)s_Q = (1 - y_Q)T_A(r_Q) + y_QT_B(r_Q)sQ=(1−yQ)TA(rQ)+yQTB(rQ)
P点的灰度级公式就为:
sP=(1−yP)TB(rP)+yPTC(rP)s_P = (1 - y_P)T_B(r_P) + y_PT_C(r_P)sP=(1−yP)TB(rP)+yPTC(rP)
问题来了: y_Q和y_P怎么获得？

从上图来看，无非就是该点到相邻的第一个子图中心的距离，所以我们需要知道第一个子图中心位置，Q点相邻的第一个子图中心坐标为(6, 2)，P点相邻的第一个子图中心坐标为(10, 2)

print(Q_numi * height_block + height_block / 2)
print(P_numi * height_block + height_block / 2)yQ = (Q[0] - (Q_numi * height_block + height_block / 2)) / height_block
yP = (P[0] - (P_numi * height_block + height_block / 2)) / height_block  # 除以height_block是为了归一化
print(yQ)
print(yP)

代码运行结果:

6.0
10.0
0.75
0.25

当然，应该也可以使用另一种方式求y_Q和y_P，但我不知道这个对不对，因为我刚开始是想着用求余的方式求y_Q和y_P，就推出了这个，和上面那种方法求出来的结果是一样的

yQ = (Q[0] - height_block / 2) % height_block / height_block
yP = (P[0] - height_block / 2) % height_block / height_block
print(yQ)
print(yP)

CLAHE代码:

def my_CLAHE(img, clipLimit, tileGridSize):height, width = img.shapeheight_block = height // tileGridSize[0]width_block = width // tileGridSize[1]total = width_block * height_blockaverage = width_block * height_block / 256ClipThreshold = int(clipLimit * average)multi_cdf = []for i in range(tileGridSize[0]):for j in range(tileGridSize[1]):grid = img[i*height_block : (i+1)*height_block, j*width_block : (j+1)*width_block]hist = calGrayHist(grid) # 统计每个子块的直方图hist = ClipHistogram(hist, ClipThreshold) # 限制对比度对直方图进行修剪# 获取累计分布直方图cdfcdf = np.zeros(256)for k in range(256):if k == 0:cdf[k] = hist[k] / totalelse:cdf[k] = cdf[k-1] + (hist[k] / total)multi_cdf.append(cdf.tolist())print(np.array(multi_cdf).shape)out = np.zeros(img.shape)# 计算CLAHE后的像素值: 根据像素点的位置，选择不同的计算方法for i in range(height):for j in range(width):# four coners，即红色区域，灰度级直接由映射函数计算得到if i <= height_block / 2 and j <= width_block / 2: # 左上角num = 0out[i][j] = int(multi_cdf[num][img[i][j]] * 255)elif i <= height_block / 2 and j >= ((tileGridSize[1] - 1) * width_block + width_block / 2): # 右上角num = tileGridSize[1] -1out[i][j] = int(multi_cdf[num][img[i][j]] * 255)elif i >= ((tileGridSize[0] - 1) * height_block + height_block / 2) and j <= width_block: # 左下角num = tileGridSize[1] * (tileGridSize[0] - 1)out[i][j] = int(multi_cdf[num][img[i][j]] * 255)elif i >= ((tileGridSize[0] - 1) * height_block + height_block / 2) and j >= ((tileGridSize[1] - 1) * width_block + width_block / 2): # 右下角num = tileGridSize[0] * tileGridSize[1] - 1out[i][j] = int(multi_cdf[num][img[i][j]] * 255)# four edges except coners，即绿色区域，灰度级由线性插值得到elif j <= width_block / 2: # 左边界num_i = int((i - height_block / 2) / height_block)num1 = num_i * tileGridSize[1]num2 = num1 + tileGridSize[1]y = (i - (num_i * height_block + height_block / 2)) / height_block # 归一化out[i][j] = int(((1 - y) * multi_cdf[num1][img[i][j]] + y * multi_cdf[num2][img[i][j]]) * 255)elif i <= height_block / 2: # 上边界num_j = int((j - width_block / 2) / width_block)num1 = num_jnum2 = num1 + 1x = (j - (num_j * width_block + width_block / 2)) / width_block # 归一化out[i][j] = int(((1 - x) * multi_cdf[num1][img[i][j]] + x * multi_cdf[num2][img[i][j]]) * 255)elif j >= ((tileGridSize[1] - 1) * width_block + width_block / 2): # 右边界num_i = int((i - height_block / 2) / height_block)num1 = num_i * tileGridSize[1] + tileGridSize[1] - 1num2 = num1 + tileGridSize[1]y = (i - (num_i * height_block + height_block / 2)) / height_block # 归一化out[i][j] = int(((1 - y) * multi_cdf[num1][img[i][j]] + y * multi_cdf[num2][img[i][j]]) * 255)elif i >= ((tileGridSize[0] - 1) * height_block + height_block / 2): # 下边界num_j = int((j - width_block / 2) / width_block)num1 = num_j + (tileGridSize[0] - 1) * tileGridSize[1]num2 = num1 + 1x = (j - (num_j * width_block + width_block / 2)) / width_block # 归一化out[i][j] = int(((1 - x) * multi_cdf[num1][img[i][j]] + x * multi_cdf[num2][img[i][j]]) * 255)# inner area，即紫色区域，灰度级由双线性插值得到else:num_i = int((i - height_block / 2) / height_block)num_j = int((j - width_block / 2) / width_block)num1 = num_i * tileGridSize[1] + num_jnum2 = num1 + 1num3 = num1 + tileGridSize[1]num4 = num2 + tileGridSize[1]x = (j - (num_j * width_block + width_block / 2)) / width_block # 归一化y = (i - (num_i * height_block + height_block / 2)) / height_block # 归一化out[i][j] = int(((1 - x) * (1 - y) * multi_cdf[num1][img[i][j]] + x * (1- y) * multi_cdf[num2][img[i][j]] + (1 - x) * y * multi_cdf[num3][img[i][j]] + x * y * multi_cdf[num4][img[i][j]]) * 255)out = out.astype(np.uint8)return out

代码运行结果:

调库实现

# 创建CLAHE对象
clahe = cv.createCLAHE(clipLimit=4.0, tileGridSize=(8, 8))
# 限制对比度的自适应阈值直方图均衡化
img_gray_clahe = clahe.apply(img_gray)

我代码中没有考虑过如果图像尺寸无法整除分块大小的情况，贴上别人的博客吧，以后有时间再去想
OpenCV自适应直方图均衡CLAHE图像和分块大小不能整除的处理

参考资源:

【图像处理算法】直方图均衡化
限制对比度自适应直方图均衡(CLAHE算法)
直方图均衡化(HE, AHE, CLAHE)
限制对比度自适应直方图均衡化（自我理解）
OpenCV–Python 图像增强（线性变换，直方图正规化，伽马变换，全局直方图均衡化，限制对比度的自适应直方图均衡化）