数学虐我千百遍
我待数学如初恋
公式打到怀疑人生
摘要
本文提出一种基于变分技术的图像感知色彩校正，提出了一个新的图像泛函，其最小值可以产生感知色彩增强后的图，这个变分公式使得局部对比度调整和数据的联系更灵活，展示了一个将梯度下降的数值实现运用到能量泛函和自动色彩增强（ACE）方程的模型。此外，欧拉-拉格朗日方程的数值近似将模型复杂度从\(O({N^2})\)减少到\(O(N\log (N))\)。
介绍
变分法适用于直方图均衡和匹配，然而最小化能量泛函进行对比度修正没有考虑到人类感知系统中的空间、局部和非线性特征，因此本文提出了关于ACE的变分公式。
本文的主要贡献：
1.考虑ACE的离散形式,可以作为某一特定能量泛函的梯度下降得到，一方面，通过研究与ACE相关的泛函推广，可知变分法可以用于基于人类感知性质的图像增强，另一方面ACE能量泛函的特殊形式揭示了离散框架中固有模型的显式特征；
2.这个公式可以更好的研究ACE的全局和局部表现，并通过不同的方式控制其影响；在变分中，ACE更容易结合数据的局部联系来避免过度增强，最后加入正则机制；
3.用梯度下降来最小化ACE泛函及欧拉-拉格朗日方程的数值近似使得模型的复杂度从\(O({N^2})\)降到\(O(N\log (N))\)。
基于变分法的全局直方图均衡
\(\psi :[0,1] \to R\)是个可微函数，考虑以下泛函：
则能量函数E(I)的变分如下：
令\(I(x) \equiv \lambda \),定义
则（2）可以写成：
若I为静态图像，且\(\delta E(I) = 0\)，则
取\(\psi’(\lambda ) = 4\lambda - 2\),若\(\psi (\lambda ) = 2{\lambda ^2} - 2\lambda + K\),K为任意常数，对于任意的\(\lambda \in [0,1]\),\({H_I}(\lambda ) = \lambda \),这就是规范直方图均衡化（uniform histogram equalization）。取K=1/2，可得\(\psi (\lambda ) = 2{(\lambda - \frac{1}{2})^2}\),代入（1）可以得到如下泛函，其最小值对应恒定的直方图均衡
如果把图像I的平均对比度记为：
因此，最小化能量函数E（I）就是最大化图像平均对比度同时最小化其偏差至理论均值1/2。
记原始图像为\({I_0}\),能量泛函的梯度下降为：
上式可求得其确定的解。使用演化等式而不是直接计算是因为（12）可以将图像\({I_0}\)的直方图逐渐转为均衡，可以随时在我们想要的时候停止，因为最后的直方图均衡结果往往看起来不真实或者过饱和。
上述式子存在两个缺点：1.演化过程保存了水平线的顺序，且没考虑到局部特征，不遵从基本的人类视觉感知特性；2.计算复杂度为\(O({N^2})\)，N为整幅图的像素个数。
图一展示了一张比较暗的图像，经过全局直方图均衡后的结果。
ACE基本形式
使用核来局部化定量，对整幅图像求和来表示，因此需要精确化图像域。将定义在\(\{ 1, \ldots W\} \times \{ 1, \ldots H\} \)的图像I扩展到\(\{ - W + 1, \ldots W\} \times \{ - H + 1, \ldots H\} \),周期为\(Z \times Z\),因此图像I的域为\({T_d}: = (Z \times Z)/(2WZ \times 2HZ)\)。
在ACE中，\(x \in {T_d}\)代表一个确定的像素target，其强度I(x)重新计算成亮度L(x)，y表示图像中不同于x的一般像素点。给定常数\(\alpha > 1\),建立如下奇函数\({S_\alpha }:[ - 1,1] \to [ - 1,1]\):
\(\alpha \)和\({S_\alpha }\)称为斜率和斜率函数
引入一个对称的权重函数\(\omega :{T_d} \times {T_d} \to {R^ + }\),且对于任意的\(x,y \in {T_d}\)都有\(\omega (x,y) = \omega (y,x)\),\(\omega (x,y)\)表示像素x,y之间的共同色彩影响权重。
用ACE计算目标像素x的亮度需要两步，第一步计算色度空间调整
然后R(x)经过动态色度再现缩放（dynamic tone reproduction scaling）变为归一化的动态范围，计算目标像素x的亮度：
其中\(M = {\max _{x \in {T_d}}}\{ R(x)\} \)
分析之前的公式发现ACE利用目标像素和其他像素之间的强度差异实现了颜色计算的空间性质。这种差异性的使用是联合实现空间对比度的可微性和灰度世界gray world(GW)原则最自然的操作。考虑一张图像只有两个像素点x和y，定义\(Dxy = I(x) - I(y)\),则\(Dxy = - Dyx\),其均值为0，由于任何奇函数的均值都为0，因此\({S_\alpha }\)的奇函数特性是ACE的重要特征。斜率函数是一个非线性变换，若强度差异较小，则经过映射可以增大对比度，反之，若强度差异较大，则函数饱和。
在ACE中，权值函数\(\omega \)实现局部性，实验证明它必须是一个单调递减函数\({\left| {x - y} \right|_{ {T_d}}}\),当\(\omega (x,y) = 1/{\left| {x - y} \right|_{ {T_d}}},y \in {T_d}\backslash \{ x\} \)时，即两个像素欧式距离的倒数，能取得较好的效果。
由于斜率函数的存在，ACE增加了色彩动态的百分比且直方图的平坦性，局部和全局对比度可以独立改变。其主要缺点就是计算复杂度为\(O({N^2})\)。ACE滤波后的图如下：
可以看到局部滤波器将暗区调整到一个更平衡且视觉效果更好的图像。
基于变分法的局部对比度增强
图像一般定义在有限域且为矩形，使用核来局部化积分表示的量，因此精确化积分域。对于给定的图像\(I:\Im \to [0,1]\),将其扩展为偶函数\({\Im ^e}: = [ - W,W] \times [ - H,H]\),在每个方向上镜面地复制图像，并以\({\Im ^e}\)为基础周期，周期性地扩展\({R^2}\)，之后图像域为\(T: = ({R^2})/(2WZ \times 2HZ)\)。
\(\omega ,\tilde \omega :T \times T \to {R^ + }\),且\(\omega (x,y) = \omega (y,x),\forall x,y \in T\)有：
\({\tilde \omega }\)也具有同样的性质。
\(J:R \to [0,\infty )\)为凸偶函数，定义图像I的平均局部对比方法如下：
平均二次局部散度如下：
其中，\({\overline I ^{\tilde \omega }}(x)\)为局部平均，有多种选择，在本文中，我们验证其中三个。
A.\({\overline I ^{\tilde \omega }}(x)\)的不同选择
1）灰度世界（Gray World）
该假设认为：对于一副有着大量色彩变化的图像，R，G，B三个分量的平均值趋于同一灰度值Gray。从物理意义上讲，灰色世界法假设自然界景物对于光线的平均反射的均值在总体上是个定值，这个定值近似地为“灰色”， 颜色平衡算法将这一假设强制应用于待处理图像，可以从图像中消除环境光的影响，获得原始场景图像。在数学上的表示如下:
上式验证了全局GW设想。GW旨在中心化图像像素值，使其位于在动态范围的中间，为了增强图像，局部地中心化信号可以最大化对比度和视觉信息。
2)加入原始数据（Attachment to Original Data）
假设\({I_0}\)为待增强的原始图像，则\({\overline I ^{\tilde \omega }}(x)\)的第二种选择如下：
定于与这种选择对应的散度为\({D_{\tilde \omega ,{I_0}}}(I)\)。在这种情况下，色散度量相当于控制图像的局部方差，在极限情况下，令\(\tilde \omega (x,y) = \delta (x - y)\),则对应的散度为
3）上述两种选择的线性结合
考虑散度是灰度世界和加入原始图像的线性结合
基于这个散度项和\({C_{\omega ,J}}\),定义泛函：
最小化能量函数(23)相当于增加对比度，同时控制所得图像的局部方差。由于最小化能量函数的目的是图像增强，因此需要增加对比度而引入噪声且不引入轮廓效果（即不过度饱和），考虑图像的视觉内容，加入原始数据可以很好地避免上述问题。
B.一般形式下能量泛函的欧拉-拉格朗日方程
对于给定输入图像\({I_0}\),我们的目标是最小化能量函数：
接下来将证明这个泛函最小值的存在性，首先计算它的变分。
使用\(\omega (x,y)\)的对称性和\(J’\)的奇函数特性，有
命题1：假设\(J\)可微，则
备注1：命题1可以拓展到任意的凸偶函数\(J\)，甚至它可以是不可微的。
\({D_{\tilde \omega ,GW}}(I)\)的初级变分为\(\delta {D_{\tilde \omega ,GW}}(I)x = 2(I(x) - (1/2))\)。对\({D_{\tilde \omega ,{I_0}}}\)有如下命题。
命题2：\({D_{\tilde \omega ,{I_0}}}(I)\)的初级变分为：
推论4.1：假设\(J\)可微，I为\({E_{\omega ,}}_{\tilde \omega }\)的最小值，则
备注2：经过适当的修改，这个推论也可以扩展到任意的凸函数\(J\)和不可微的情况。
C.\({E_{\omega ,}}_{\tilde \omega }\)最小值的存在性
实际上，\({E_{\omega ,}}_{\tilde \omega }\)最小值的存在性可以通过一些附加的假设和离散形式进行证明。将能量函数写成离散形式，假设\[I,{I_0}:{T_d} \to {R^ + }\]表示两幅图像，其中\({I_0}\)为原始图像，\(\omega ,\tilde \omega :{T_d} \times {T_d} \to {R^ + }\)为对称函数，则
考虑离散泛函：
命题3：若函数\(J\)对于任意的常数C都满足\(J(r) \le C|r\)，则能量函数\({E^d}_{\tilde \omega ,\omega }\)对于满足条件的图像\(I:{T_d} \to {R^ + }\)有最小值。这个命题证明了\({E^d}_{\tilde \omega ,\omega }\)最小值的存在性。
备注3：若用\({C_{\omega ,J}}(G * I)\)替代\({C_{\omega ,J}}(I)\)，在连续域也有相似的命题。也就是说，用正则化后的形式\(G * I\)代替图像I，其中G为卷积核，如高斯。
ACE变分公式
本部分主要介绍在上一章节中能增强局部对比度的变分技术也可以得到ACE局部色彩校正均衡。为了与之前的离散框架建立连接，使用离散能量函数。
其中，A(x)为衡量A，即对任意的\(x \in {T_d}\)都有\(A(x) = A\)，结合归一化因子(14)来定义权重，如下：
有了上述权重，(14)对应的ACE空间色度对比可以些成：
因此
目标是使用上一章节中采用的变分法将I(x)转化为L(x)。
考虑如下泛函：
其中，\(S{‘_\alpha } = {s_\alpha }\)，\(S{‘_\alpha }\)为凸偶函数，其导数为奇函数，\({s_\alpha }\)为ACE斜率函数且\(0 < M \le 1\)是之后确定的一个常数，每一项前的系数为正则化常数。
若I为E的最小值，那么\(\delta E(I) = 0\)，结合推论4.1可知
等式右侧为ACE计算得到的亮度L(x)，因此上式表明L(x)为ACE的一个确定点。
备注4：若采用梯度下降策略来最小化能量函数\(E(I)\)，需要解决的是：
若使用显格式（explicit scheme）离散化参数t，可得
一方面，令\(M: = {\max _{x \in {T_d}}}\{ {R_{ {I^0}}}(x)\} \),\(\Delta t = 1,k = 0\)可得：
另一方面，用(42)作为数值法得到稳定状态的ACE，但这个稳定状态往往会被过度增强干扰，使得输出结果视觉效果差。因此使用散度方式(22)，但其会偏离原始图像，需要附加数据项。为了简化，在实验部分，使用\(\tilde \omega (x,y) = \delta (x - y)\)，因此\({\overline I ^{\tilde \omega }}(x) = \sum\nolimits_{z \in {T_d}} {\tilde \omega (x,y),} {I^{_0}}(y) = {I^{_0}}(x)\),参数\(\beta = (1/2),\gamma = \lambda /2,\lambda > 0\)，在这种情况下，梯度下降方程为：
化简后得
其中，\(0 < \Delta t \le 1/(1 + \lambda )\),等式右侧包含一项在相反方向上进行平衡来抑制与原始图像的偏离。
A.均匀直方图均衡与变分框架中ACE的比较
对比均匀直方图均衡和ACE变分公式，两个泛函分别为
两个泛函中第一项都包含一个中间灰度值附近的图像信号能量分布，最小化这部分对整幅图产生局部影响。第二项不同：对于均匀直方图均衡，利用两个像素间的绝对差异可得全局对比度方法，但不适合感知视角；ACE变分公式中，采用了更复杂的对比度方法，由于斜率和权值函数的引入，可以重建感知效果，其非线性和局部特征更适合人类的感知系统（HSV）。
B.模型的重要性质：局部感知和水平线
根据形态学，图像I的特征在于其上（或下）水平集\({X_\lambda } = \{ x:I(x) \ge \lambda \} (resp.X{‘_\lambda } = \{ x:I(x) \le \lambda \} )\)对于任意的\(\lambda \)上（下）水平集由图像像素值在上(下)面的点组成，\(\lambda \in R\)。图像I可以用下面的重建公式从水平集中重建得到：
水平线通常被定义为水平集的边界，并且直观地表示图像的等级线。需要研究的是水平集对非线性局部色彩增强算法的影响。从(45)中分析这一操作的性质:
首先分析全局对比度增强的情况。
命题4:\(J:R \to [0,\infty )\)为凸函数，\(\omega (x,y) = 1/Area(T)\),\({H_J} = {H_{\omega ,J}}\),则
在标量情况下，若\(I(x) < I(\tilde x)\),则\({H_J}(I(x)) \le {H_J}(I(\tilde x))\)，若\(J\)为严格凸函数则不等式也严格成立。根据(45)和\(\Delta t < 1/(1 + \lambda )\)，若\({I^k}(x) < {I^k}(\tilde x)\)，则\({I^{k + 1}}(x) < {I^{k + 1}}(\tilde x)\)。另一方面，若\(I(x) = I(\tilde x)\),则\({H_J}(I(x)) = {H_J}(I(\tilde x))\)。因此若\({I^k}(x) = {I^k}(\tilde x)\)，则\({I^{k + 1}}(x) = {I^{k + 1}}(\tilde x)\)。这表明从\({I^k}\)到\({I^{k + 1}}\)的映射可以转变为从\({I^k}\)到\({I^{k + 1}}\)的水平集，在全局情况下（\(\omega = 1/Area(T)\)）算法没有破坏\({I_0}\)的映射结构。换句话说，对于任何迭代k，\({I^k}\)的水平线族与\({I_0}\)的水平线族一致。若包含卷积核\(\omega (x,y)\)就违背了上述性质，但只影响图像I中一小部分平滑区域，对局部影响不明显。
总的来说，这一算法保持了图像的局部关系，同时全局上有所改变。这个公式满足基本的视觉感知特性（局部性和水平线的顺序）。
数值近似（优化）
由于斜率函数\({s_\alpha }\)不可微，为了更好地近似，用sigmoid函数进行替换，记为\({ {\tilde s}_\alpha }\)。考虑离散图像\(I:{T_d} \to [0,1]\)为\({R^{4WH}}\)的向量，定义
可以看到，对任意的\(x \in {T_d}\)，都有\(|{ {\tilde R}_I}(x)| \le 1\)。便于估计，考虑M=1为\({ {\tilde R}_I}(x)\)的上界，则数值法(45)可写成
因为\({I^0}\)在[0,1]，若\({I^k}\)也属于[0,1],则\(|{ {\tilde R}_{ {I^k}}}(x)| \le 1\)且
因此对于任意的迭代\(k \ge 0\)，\({I^k}\)都属于[0,1]。由于在离散形式下，这表明在\({\{ {I^k}\} _k}\)中存在可收敛的子序列，这是一种弱形式的稳定性。 另一方面，对于一些常数C>0,我们有
为了证明上式，定义\(F(t) = E({I^k} + t({I^{k + 1}} - {I^k}))\)，使用中值定理，存在\({t^{*} } \in (0,1)\)使得\(F(1) - F(0) = F’({t^{*} })\)。根据(52)有\(\left| { {I^{k + 1}} - {I^k}} \right| \le 2(1 + \lambda )\Delta t\left| {\rm{1}} \right|\)（其中，1表示连续图像等于1）因此:
其中，\({C_{\rm{0}}}{\rm{ > 0}}\)为E二阶导的范围（依赖于\({ {\tilde s}_\alpha }\)二阶导的范围）。式子(53)表明
当\(\Delta t \to 0\)时，右边趋于0，意味着能量函数的导数为负。可以得出结论，我们的数值方案是稳定的（在某种意义上说它不会发散）并将能量降低到有序误差\({(\Delta t)^2}\)。 虽然不能证明(52)完全收敛，但在所有的实验中，不到60次迭代中达到了稳定状态。
备注5:使用一个确定点论证能很容易地证明当\(\lambda > {\max _{r \in [ - 1,1]}}|(d{ {\tilde s}_\alpha }/dr)(r)| - 1\)时，数值法(52)的收敛性。这样表明与数据的联系更强了。
A.算法加速
根据上式可知，ACE的计算复杂度为\(O(N)\)乘以N个像素，所以总的复杂度为\(O({N^2})\)。因此使用数值近似来加速，第一步，将函数\({ {\tilde s}_\alpha }\)近似为一个多项式(在有限域中)可以将复杂度减少到确定数量的卷积；第二步，使用快速傅里叶变换(FFT)将最终的计算复杂度减少到\(O(N\log (N))\)。
第一步，将函数\({ {\tilde s}_\alpha }\)近似为一个[-1,1]的多项式，如下
即最小化下式：
其中，\(p(r)\)为一类度为n的多项式。
使用多项式来估计\({ {\tilde s}_\alpha }^{(n)}(I(x) - I(y))\)
观察上式，发现中括号内的部分只与\(I(x)\)和系数\({c_i}\)有关，把它记为\({a_i}\),则\({ {\tilde s}_\alpha }^{(n)}(I(x) - I(y)) = \sum\limits_{i = 0}^n { {a_i}} I{(y)^i}\)，因此可得\(R(x)\)的近似结果如下：
前一项的计算复杂度为\(O(1)\),后一项使用FFT，计算复杂度为\(O(N\log (N))\)
算法调整
- 改变斜率函数调整最后的对比度：斜率越大，对比度越大。
- 与基本ACE一样，核的形状及其宽度是造成像素间相互色彩影响的重要因素。滤波器局部形状的直观例子如图5，选择核的形状是一种控制图像校正局部影响的方式。
- 在ACE的最终形式中，加入两个参数Keep Original Grey (KOG) and Keep Original Dynamic
Range (KODR)。KOG保留原始图像均值，KODR保持其原始动态范围，防止图像过拟合，更接近原始图像。实验
参数设置：\(\Delta t{\rm{ = 0}}{\rm{.15,}}\lambda {\rm{ = 1,n = 7,}}\alpha {\rm{ = 10}}\),\(\omega \)为高斯核。
图4展示了一张灰度图像经过全局均衡和我们的算法后的结果，可以发现本文提出的算法根据其到边缘的距离对灰度值有较好的局部修正，更符合HVS感知特性。
图6显示了该算法在图像去量化中的应用。 左上角的原始图像，编码为8位，256色。 左下方红色通道直方图显示典型的不相交尖峰，而右侧重新计算中间值显示，贯穿整个24位调色板。
总结
本文我们主要研究了基于变分法的全局和局部色彩校正算法，尤其是考虑到人类色彩感知重要特征的色彩均衡模型ACE。此外，使用多项式近似将算法复杂度从\(O({N^2})\)减少到\(O(N\log (N))\)。
