OpenCV—python 自动色彩均衡（ACE）

文章目录

一、ACE算法理论
- 1.1 ACE算法
二、快速ACE算法
- 2.1 原理
- 2.2 代码演示

一、ACE算法理论

图像对比度增强的算法在很多场合都有用处，特别是在医学图像中，这是因为在众多疾病的诊断中，医学图像的视觉检查时很有必要的。医学图像由于本身及成像条件的限制，图像的对比度很低。因此，在这个方面已经开展了很多的研究。这种增强算法一般都遵循一定的视觉原则。众所周知，人眼对高频信号（边缘处等）比较敏感。虽然细节信息往往是高频信号，但是他们时常嵌入在大量的低频背景信号中，从而使得其视觉可见性降低。因此适当的提高高频部分能够提高视觉效果并有利于诊断。

ACE在图像处理方面可以有两种，
一种是：Automatic Color Equalization，即自动彩色均衡。【论文与代码】
一种是：Adaptive Contrast Enhancement，即自适应对比度增强

高动态范围图像是指在一幅图像中，既有明亮的区域又有阴影区域，为了使细节清晰，需要满足以下几点：

对动态范围具有一定的压缩能力；
对亮暗区域的细节有一定的显示能力；
满足上面条件基础上不破坏图像的清晰度。

对于高动态范围处理，基于人眼视觉系统（HSV）在颜色连续和亮度连续方面得到较好的满足。该算法考虑了图像中颜色和亮度的空间位置关系，进行局部的自适应滤波，实现具有局部和非线性特征的图像亮度，色彩与对比度调整，同时满足灰度世界理论和白斑点假设。

算法原理
Rizzi等依据Retinex理论提出了自动颜色均衡算法，该算法考虑了图像中颜色和亮度的空间位置关系，进行局部特性的自适应滤波，实现具有局部和非线性特征的图像亮度与色彩调整和对比度调整，同时满足灰色世界理论假设和白色斑点假设。

1.1 ACE算法

获得空域重构图像
对图像进行色彩/空域调整，完成图像的色差矫正，得到空域重构图像：

Rc=∑j∈subset,j≠pr(Ic(p)−Ic(j))d(p,j)(1)R_c = \sum_{j\in subset,j \neq p }\frac{r(I_c(p) - I_c(j))}{d(p,j)} \tag{1}Rc=j∈subset,j=p∑d(p,j)r(Ic(p)−Ic(j))(1)

其中Rc(p)Rc(p)Rc(p)是中间结果，Ic(p)−Ic(j)I_c(p) - I_c(j)Ic(p)−Ic(j)为2个点的亮度差，d(p,j)d(p,j)d(p,j) 表示距离度量函数，r(∗)r(∗)r(∗) 为亮度表现函数，需要是奇函数，这一步可以适应局部图像对比度，r(∗)r(∗)r(∗) 可以放大较小的差异，并丰富大的差异，根据局部内容扩展或者压缩动态范围。一般的，r(∗)r(∗)r(∗) 为:
r(x)={1x<−Tx/T−T⩽x⩽T−1x>Tr(x) = \left\{\begin{matrix} 1 & x<-T\\ x/T & -T \leqslant x \leqslant T\\ -1 & x>T \end{matrix}\right.r(x)=⎩⎨⎧1x/T−1x<−T−T⩽x⩽Tx>T
对校正后的图像进行动态扩展。ACE算法是对单一色道进行的，对于彩色图片需要对每一个色道分别处理。一种简单的线性扩展可以表示为：
Oc(p)=round[127.5+scRc(p)](2)O_c(p) = round[127.5+s_cR_c(p)] \tag{2}Oc(p)=round[127.5+scRc(p)](2)
其中：scs_csc斜率为：[(mc,0),(mc,255)][(m_c,0),(m_c,255)][(mc,0),(mc,255)]，Mc=max⁡p[Rc(p)]，mc=min⁡p[Rc(p)]M_c = \max_{p}[R_c(p)]，m_c = \min_{p}[R_c(p)]Mc=maxp[Rc(p)]，mc=minp[Rc(p)]

还可以将其映射到 [0,255][0,255][0,255]的空间中：
L(x)=R(x)−min⁡Rmax⁡R−min⁡RL(x) = \frac{R(x) - \min R}{ \max R - \min R}L(x)=maxR−minRR(x)−minR

自动彩色均衡算法改进:
式（1）算法复杂度较高，对于一副像素数为N的图像，需要执行O(N2)级次非线性映射计算，图像尺寸越大，耗时越多，所以针对式（1）产生了许多加速改进算法。例如：LLLUT加速策略，使用快速傅里叶变换替换卷积，将ACE转换为对规范直方图均衡化的一种平滑和局部修正的方法，并给出了求解最优模型：
argmin⁡I12∑x(I(x)−12)2−14M∑x∑y≠xω(x,y)Sα(I(x)−I(y))arg \min_I\frac{1}{2}\sum_{x}(I(x)-\frac{1}{2})^2 -\frac{1}{4M}\sum_x\sum_{y\neq x}\omega(x,y)S\alpha (I(x)-I(y))argImin21x∑(I(x)−21)2−4M1x∑y=x∑ω(x,y)Sα(I(x)−I(y))

其中：: Sα′=sαS_{\alpha }' = s_\alphaSα′=sα
ω=1/∣∣x−y∣∣\omega =1/||x-y||ω=1/∣∣x−y∣∣
M=max⁡xR(x)M = \max_xR(x)M=maxxR(x)

对于改进方法，可以考虑的因素：

其他的坡度函数SαS \alphaSα，多项式函数逼近；
除了 1/∣∣x−y∣∣1/||x-y||1/∣∣x−y∣∣ 外的权重函数的选择；
在求和的过程中，y 可以限制在一个小窗口中；
L(x)L(x)L(x) 的一些其他的方法；

多项式逼近展示：

二、快速ACE算法

2.1 原理

ACE的增强效果普遍与retinex好。需要注意的是，ACE中当前像素是与整个图像的其他像素做差分比较，计算复杂度非常非常高，这也是限制它应用的最主要原因。
本文主要基于两个假设：
（1）对一副图像ACE增强后得到输出Y，如果对Y再进行一次ACE增强，输出仍然是Y本身；
（2）对一副图像的ACE增强结果进行尺寸缩放得到Y，对Y进行ACE增强，输出仍然是Y本身。

如果上面假设成立，我们就可以对图像进行缩放得到 I1I_1I1，对 I1I_1I1 的ACE增强结果进行尺度放大（与I尺寸一样）得到Y1Y_1Y1，那么 YYY 和 Y1Y_1Y1 是非常接近的，我们只需要在Y1基础上进一步处理即可。
这里就又引申了两个细节问题：

如何快速的求 I1I_1I1 的ACE增强结果？
其实很简单，对它再次缩放得到 I2I_2I2 ，求 I2I_2I2 的增强结果，依次类推，就是金字塔结构思想。
如何在 Y1Y_1Y1 基础上进一步处理得到 YYY？
因为是在整个图像域进行差分比较运算，与近处邻域像素的比较构成了YYY的细节信息，与远处像素的比较构成了YYY的全局背景信息，那么我们合理假设，YYY和Y1Y_1Y1的全局背景信息相同，只更新细节信息即可，也就是，我们需要在Y1Y_1Y1基础上加上 III 中邻近像素的差分结果，并减去Y1Y_1Y1中邻近像素的差分结果就是最终的输出 YYY。

为叙述方便，这里假设后面的图像都是归一化到 [0,1] 之间的浮点数图像。
ACE算法的计算公式为：
Y=∑g(I(x0)−I(x))∑w(x0,x)Y = \frac{ \sum g(I(x_0)-I(x))}{\sum w(x_0,x)}Y=∑w(x0,x)∑g(I(x0)−I(x))

其中：: x∈Ix\in Ix∈I
www是权重参数，离中心点像素越远w值越小，可以直接取值欧氏距离。
g()g()g()是相对对比度调节参数，非线性的，简单取如下计算方法：
g(x)=max⁡(min⁡(ax,1.0),−1.0)g(x) = \max(\min(ax, 1.0), -1.0)g(x)=max(min(ax,1.0),−1.0)
a是控制参数，值越大，细节增强越明显。计算完后，还要对Y进行一次归一化即可得到最终的增强图像。

2.2 代码演示

import cv2
import numpy as np
import mathdef stretchImage(data, s=0.005, bins = 2000):    #线性拉伸，去掉最大最小0.5%的像素值，然后线性拉伸至[0,1]ht = np.histogram(data, bins);d = np.cumsum(ht[0])/float(data.size)lmin = 0; lmax=bins-1while lmin<bins:if d[lmin]>=s:breaklmin+=1while lmax>=0:if d[lmax]<=1-s:breaklmax-=1return np.clip((data-ht[1][lmin])/(ht[1][lmax]-ht[1][lmin]), 0,1)g_para = {}
def getPara(radius = 5):                        #根据半径计算权重参数矩阵global g_param = g_para.get(radius, None)if m is not None:return msize = radius*2+1m = np.zeros((size, size))for h in range(-radius, radius+1):for w in range(-radius, radius+1):if h==0 and w==0:continuem[radius+h, radius+w] = 1.0/math.sqrt(h**2+w**2)m /= m.sum()g_para[radius] = mreturn mdef zmIce(I, ratio=4, radius=300):                     #常规的ACE实现para = getPara(radius)height,width = I.shapezh,zw = [0]*radius + range(height) + [height-1]*radius, [0]*radius + range(width)  + [width -1]*radiusZ = I[np.ix_(zh, zw)]res = np.zeros(I.shape)for h in range(radius*2+1):for w in range(radius*2+1):if para[h][w] == 0:continueres += (para[h][w] * np.clip((I-Z[h:h+height, w:w+width])*ratio, -1, 1))return resdef zmIceFast(I, ratio, radius):                #单通道ACE快速增强实现height, width = I.shape[:2]if min(height, width) <=2:return np.zeros(I.shape)+0.5Rs = cv2.resize(I, ((width+1)/2, (height+1)/2))Rf = zmIceFast(Rs, ratio, radius)             #递归调用Rf = cv2.resize(Rf, (width, height))Rs = cv2.resize(Rs, (width, height))return Rf+zmIce(I,ratio, radius)-zmIce(Rs,ratio,radius)    def zmIceColor(I, ratio=4, radius=3):               #rgb三通道分别增强，ratio是对比度增强因子，radius是卷积模板半径res = np.zeros(I.shape)for k in range(3):res[:,:,k] = stretchImage(zmIceFast(I[:,:,k], ratio, radius))return resif __name__ == '__main__':m = zmIceColor(cv2.imread('p4.bmp')/255.0)*255cv2.imwrite('zmIce.jpg', m)

鸣谢
文献地址：https://www.researchgate.net/publication/253622155_Real-Time_Adaptive_Contrast_Enhancement_For_Imaging_Sensors
C++:https://blog.csdn.net/weixin_45709330/article/details/104447450
论文与源码：http://www.ipol.im/pub/art/2012/g-ace/
https://blog.csdn.net/zmshy2128/article/details/53470357#comments
https://www.cnblogs.com/Imageshop
https://blog.csdn.net/weixin_45709330/article/details/104447336