总变差(Total Variation, TV)正则[1]在变分法图像处理中具有广泛的应用, 其变分形式为如下的最优化问题:
$ \mathop {\min }\limits_u \int_\mathit{\Omega } {\left| {\nabla u(x)} \right|{\rm{d}}x} $ | (1) |
它通过惩罚导数来保持图像边缘的不连续性, 具体来说就是图像能量变分模型在随人工时间的演化过程中, 图像能量只沿着图像梯度正交的方向进行扩散。这种能量扩散的各向异性的优点是在有效保护图像边缘强度的同时沿图像边缘进行平滑从而避免重建图像中的锯齿现象,但是其缺点在于TV正则问题解的分段常数性质容易造成在重建图像平坦区域产生分块现象[2]。
把带有噪声的图像看作热传导方程ut=Δu在零时刻的初始条件, 以这个初值问题的解作为图像滤波的结果, 能够达到很好的去噪结果[3]。热传导方程是下面的最小化问题的稳定状态解:
$ \mathop {\min }\limits_u \frac{1}{2}\int_\mathit{\Omega } {{{\left| {\nabla u(x)} \right|}^2}{\rm{d}}x} $ | (2) |
这是一种各向同性的扩散, 能够平滑噪声、消除分块现象,但是随着人工时间的增大, 在去噪效果越好的同时边缘也越来越模糊。图像边缘被认为是图像的最重要的特征, 而且也是人的视觉最为敏感的特征, 所以这是一个让人无法接受的结果。
为了既平滑图像平坦区域又保护图像边缘, 一个自然的想法是把这两个方程合二为一, 自适应地在式(1) 和式(2) 之间转换, 即在图像平坦区域式(2) 其作用, 在图像边缘式(1) 起作用。基于此, Blomgren等[4]提出了下面的变指数正则模型:
$ \mathop {\min }\limits_u \int_\mathit{\Omega } {{{\left| {\nabla u(x)} \right|}^{p(\left| {\nabla u} \right|)}}{\rm{d}}x} $ | (3) |
Chen等[5]用分段函数定义变指数函数p, 并研究了这个模型的解的存在性、唯一性以及扩散性质。Huang等[6]在变指数p-Laplace能量泛函中融入区域信息, 提出一个活动轮廓模型用于图像分割。在这个模型中变指数能量泛函正则零水平集, 使水平集向精确的目标边界运动。Liu等[7]研究了带有变指数Laplace的反应扩散方程, 并用于图像去噪。Maiseli等[8]获得了变指数扩散驱动正则泛函用于多帧图像超分辨率重建。它能根据图像局部特征自适应的调整扩散机制, 在图像平坦区域进行线性各向同性扩散而在图像边缘或轮廓扩散消失。董婵婵等[9]提出了一种基于变指数的片相似性扩散图像降噪算法。它在变指数自适应降噪模型基础上, 引入片相似性, 获得了一个新的边缘检测算子和扩散系数函数。张芳等[10]提出一个变指数和非局部的最大似然期望最大低剂量CT重建算法。该算法考虑到传统各向异性扩散中降噪的不充分性, 有效融合热传导和各向异性扩散到变指数函数中, 把它和代替梯度检测边缘的相似度函数运用到传统各向异性扩散中, 从而达到所期望的效果。王益艳[11]利用图像局部特征提出了一种基于Lp范数的变指数正则变分模型。这个模型采用结构张量作为Lp范数算子的自适应调整参数, 克服传统算子对噪声敏感的缺陷。
本文通过减小插值图像边缘宽度来抑制插值图像模糊, 通过分析减小图像边缘宽度与图像能量扩散之间的关系, 给出了一个新的变指数函数。它能在靠近图像边缘中心附近进行各向异性的能量扩散, 保持图像边缘的清晰度, 而在离图像边缘中心较远的区域进行各向同性的能量扩散, 消除插值边缘的锯齿现象并消除图像平坦区域的分块现象。数值实验结果表明无论是从主观的视觉评价, 还是从客观的全局性能评价(平均结构相似度), 本文方法既能很好地重建插值图像的边缘, 又不会在插值图像中产生锯齿现象以及分块现象。
1 p-Dirichlet泛函与变指数泛函广泛应用于图像处理中的p-Dirichlet泛函具有如下形式[12]:
$ E(u) = \int_\mathit{\Omega } {\frac{1}{p}{{\left| {\nabla u} \right|}^p}{\rm{d}}\mathit{\Omega }} $ | (4) |
Ω是Rn中的具有Lipshitz边界的有界开子集。在适当的边界条件下, 最小化E可以得到相应的p-Laplacian方程δE=0:
$ \delta E =-\nabla ({\left| {\nabla u} \right|^{p-2}}\nabla u) $ | (5) |
上述p-Laplacian方程的稳定状态解用最速下降法可表示为:
$ {u_t} = \nabla ({\left| {\nabla u} \right|^{p-2}}\nabla u) $ | (6) |
对p=1来说, 式(4) 变为经典的TV变分模型, 对应的Euler-Lagrange方程为:
$ {u_t} = \kappa \left| {\nabla u} \right| $ | (7) |
其中κ是局部梯度意义下的欧拉曲率。式(7) 的稳定状态解的迭代求解过程可以看作在人工时间t下的能量扩散。这个扩散具有各向异性扩散性质, 它使图像能量沿着图像梯度正交方向扩散, 避免了对图像边缘的模糊。从物理运动的观点可以解释为, 图像能量在式(7) 的作用下, 在t的迭代过程中会保持图像轮廓的位置和强度从而保持图像边缘的清晰度。
对p=2来说, 式(6) 变为热传导方程ut=Δu。上述方程的能量演化具有各向同性的性质, 能够平滑掉图像噪声, 但是它也使得图像边缘变得模糊, 不适合图像重建。p→ ∞的情况就是众所周知的无穷Laplace。
在图像去噪、重建等领域, 往往需要在图像平坦区域进行平滑以去噪同时保持图像边缘的清晰度。TV变分模型以及热传导方程都只能起到单一的作用, 兼顾二者的算子是现实需要的。Blomgren等[4]提出了下面的最小化问题:
$ \mathop {\min }\limits_{u \in BV \cap L(\mathit{\Omega })} \int_\mathit{\Omega } {{{\left| {\nabla u} \right|}^{p(\left| {\nabla u} \right|)}}{\rm{d}}x} $ | (8) |
其中
要想通过式(8) 获得满意的图像处理效果, 适当地选择单调下降函数p是关键。下面通过对图像边缘的几何分析探索函数p应具有的性质。如果把数字图像的边缘看作是类斜面的剖面, 那么图像边缘越模糊则斜坡坡度越小, 对应的图像边缘宽度越大;反之则斜坡坡度越大而图像边缘宽度越小。因此, 要使重建图像有清晰的边缘就需要减小图像边缘的宽度使斜坡坡度更大。增加斜坡坡度可以通过降低斜坡下半部灰度值同时增加斜坡上半部灰度值, 如图 1所示。把图像像素灰度值看作质量粒子, 斜坡上像素灰度值的变化可以通过图像能量扩散推动质量粒子运动而实现。可以推动斜坡顶上的质量粒子加速向斜坡中心方向扩散, 使质量粒子大量堆积在斜坡上半部;同时推动斜坡中心下半部质量粒子加速向斜坡底扩散, 使质量粒子分散在斜坡坡底。这样在图像边缘斜坡中心较窄的宽度范围内灰度对比度变大, 实现减小图像边缘宽度的目的, 从而形成清晰的图像边缘。这就要增强在图像边缘两边的区域的各向同性扩散(p≈2), 而且增强在靠近边缘的区域的各向同性扩散(p≈1)。由此, 函数p可以选取如下形式:
$ p(x) = 1 + \exp (-{(\varepsilon x)^\alpha }) $ | (9) |
图 2为函数p取不同参数时的图像。可以看出:x比较小的时候, 函数p趋近于2;x比较大的时候, 函数p趋近于1。函数p的两个参数α, ε起着对图像像素分类的作用, 即多大梯度的像素需要执行各向同性的扩散(平滑), 相应的其他像素需要各向异性的扩散。鉴于以上考虑, 本文提出如下的图像插值能量泛函:
$ \mathop {\min }\limits_{u \in BV \cap L(\mathit{\Omega })} \int_\mathit{\Omega } {\varphi (x, \left| {\nabla u} \right|){\rm{d}}x} + \frac{\lambda }{2}{\left| {Hu-{u_0}} \right|^2} $ | (10) |
$ \varphi (x, \left| {\nabla u} \right|) = {\left| {\nabla u} \right|^{1 + \exp (-{{(\varepsilon \left| {\nabla u} \right|)}^\alpha })}} $ | (11) |
λ是正实数。
在这个模型中图像中心区域的远端图像梯度较小, 此时函数p趋近于2, 达到增强各向同性扩散的能力;反之, 在图像中心区域的近端图像梯度较大, 此时函数p趋近于1, 达到增强各向异性扩散的能力。这样就能减小图像边缘宽度, 增强图像边缘两边像素灰度的比率, 从而产生清晰的图像边缘。
这个能量泛函的Euler-Lagrange方程为:
$ \frac{{\partial u}}{{\partial t}}- {\rm{div}}(\varphi (x, |\nabla u|)) + \boldsymbol{\lambda} {\boldsymbol{H}^{\bf{\boldsymbol{T}}}}(\boldsymbol{H}u- {u_0}) = 0;{\rm{ }}\mathit{\Omega } \times [0, T] $ | (12) |
$ \frac{{\partial u}}{{\partial n}}(x, t) = 0;{\rm{ }}\partial \mathit{\Omega } \times [0, T] $ | (13) |
$ u(0) = 0;\mathit{\Omega } $ | (14) |
其中:矩阵H用来模拟图像获取过程中滤波和下采样过程,u是原始的高分辨率图像。
本文采用显式的有限差分格式对式(12) 是进行迭代计算, 其扩散项:
$ \begin{array}{l} {\rm{div}}(\varphi ) = |\nabla u{|^{p(|\nabla u|)- 2}}[(p(|\nabla u|)-1)\Delta u + \\ (2-p(|\nabla u|))|\nabla u|\kappa + \nabla p \cdot \nabla u\log |\nabla u|] \end{array} $ | (15) |
本章用式(15) 对自然图像、纹理图像进行插值, 以说明提出算法的有效性。插值结果显示在图 3~5中。插值输入图像是原始图像用Matlab函数imresize获得, 它包含了经过低通滤波和下采样过程。这些输入图像用上一章的算法再恢复到原始图像尺寸大小从而实现插值。同时, 以Chen等[5]提供的方法以及运用加权最小二乘法的鲁棒软决策插值(Robust Soft-decision Adaptive Interpolation, BSAI)方法[13]的实验结果作为对比实验。BSAI方法的实验结果由作者提供的软件实现。在实验中本文算法中的α=4, 时间步长Δt=0.15, 其他实验参数的优选基于观察者对图像边缘是否清晰、是否有锯齿现象、在平坦区域或边缘附近是否有振铃现象等的主观评价作出。
与其他算法相比较, 本文算法真实地重构了低分辨率图像中的细微信息, 这正是图像插值问题所希望获得的效果。
图 3中, 本文算法重构的背景纹理(原始图像最左边区域)看起来比其他算法获得的结果更正确、更清晰, 纹理也更丰富。特别是在原始图像中间箭头指示的阴影区域的纹理, Chen和BSAI算法的重建结果几乎看不出有纹理, 而本文算法结果中纹理依然可辨。在对帽檐的层层线条状纹理(原始图像上边的箭头指示的条纹处)的重建中, 本文算法结果也是最清晰的, BSAI算法结果清晰度最弱。帽子下面的木板纹理在本文算法结果中清晰度也是最高的, 特别是原始图像最右边箭头指示区域。
在图 4中BSAI算法和本文算法结果中赛车手较清晰一些, 文献[5]结果稍微有些模糊。赛车轮胎齿可以清晰反映出本文算法在保持图像边缘方面的优势, 本文算法产生了更清楚的轮胎齿, BSAI结果可辨性较差, 而文献[5]结果几乎不可辨。赛车轮胎的辐条也具有同样的可辨性。
图 5中, 雕塑图像人物头像左边鬓角处的头发, 本文算法结果看起来更自然, 更真实;而BSAI算法结果较模糊, 文献[5]算法结果在这个地方几乎不会引起视觉上的注意。从图 5中的球也可以看出, 本文算法对球面上的纹路保持得最好, 其他两种算法对原始图像的纹路平滑得更严重一些。
平均结构相似度指标(Mean Structural SIMilarity, MSSIM)[14]在这里用来刻画原始参考图像与插值输出图像之间的差异。MSSIM指标比峰值信噪比或其他指标更能刻画图像的视觉质量的好坏。MSSIM指标取值为[0, 1]区间, 值越大图像的视觉质量越好, 其计算原始代码参见http://www.cns.nyu.edu/$\sim$lcv/ssim/(用缺省的参数)。本文所使用的测试图像有帽子图像、赛车图像、雕塑图像、花以及船图像, 通过计算MSSIM, 从客观指标上对3种算法进行比较。表 1显示了几种算法对测试图像计算的MSSIM, 从表 1中可以看出本文方法的客观指标在所有实验中都有明显的改善。这表明本文模型具有较好的效果。
本文提出了一种结合TV变分和热扩散的变指数变分模型图像插值方法。这种模型通过分析图像插值变分模型的扩散特性, 定义了一个指数函数。这个指数函数使图像能量在随人工时间的演化过程中在图像边缘附近只沿着图像轮廓扩散, 以消除图像边缘在插值过程中的震荡进而获得光滑的边缘;在图像平坦区域进行热扩散消除分块现象等人工虚像。
[1] | RUDIN L I, OSHER S, FATEMI E. Nonlinear total variation based noise removal algorithms[J]. Physica D:Nonlinear Phenomena, 1992, 60(1/2/3/4): 259-268. |
[2] | RING W. Structural properties of solutions to total variation regularization problems[J]. Mathematical Modelling and Numerical Analysis, 2000, 34(4): 799-810. DOI:10.1051/m2an:2000104 |
[3] | WEICKET J. Anisotropic Diffusion in Image Processing[M]. Stuttgart: Teubner-Verlog, 1998: 59-60. |
[4] | BLOMGREN P, CHAN T F, MULET P, et al. Total variation image restoration:numerical methods and extensions[C]//Proceedings of the 1997 International Conference on Image Processing. Piscataway, NJ:IEEE, 1997:384-387. |
[5] | CHEN Y M, LEVINE S, RAO M. Variable exponent, linear growth functionals in image restoration[J]. SIAM Journal on Applied Mathematics, 2006, 66(4): 1383-1406. DOI:10.1137/050624522 |
[6] | HUANG C, ZENG L. Level set evolution model for image segmentation based on variable exponent p-Laplace equation[J]. Applied Mathematical Modelling, 2016, 40(17/18): 7739-7750. |
[7] | LIU Q, GUO Z, WANG C. Renormalized solutions to a reaction-diffusion system applied to image denoising[J]. Discrete and Continuous Dynamical Systems:Series B, 2016, 21(6): 1829-1858. |
[8] | MAISELI B J, ELISHA O A, GAO H. A multi-frame super-resolution method based on the variable exponent nonlinear diffusion regularizer[J]. EURASIP Journal on Image and Video Processing, 2015, 2015(1): 22. DOI:10.1186/s13640-015-0077-2 |
[9] | 董婵婵, 张权, 郝慧艳, 等. 基于变指数的片相似性扩散图像降噪算法[J]. 计算机应用, 2014, 34(10): 2963-2966. (DONG C C, ZHANG Q, HAO H Y, et al. Patch similarity anisotropic diffusion algorithm based on variable exponent for image denoising[J]. Journal of Computer Applications, 2014, 34(10): 2963-2966. DOI:10.11772/j.issn.1001-9081.2014.10.2963) |
[10] | 张芳, 崔学英, 张权, 等. 基于变指数各向异性扩散和非局部的最大似然期望最大低剂量CT重建算法[J]. 计算机应用, 2014, 34(12): 3605-3608. (ZHANG F, CUI X Y, ZHANG Q, et al. MLEM low-dose CT reconstruction algorithm based on variable exponent anisotropic diffusion and non-locality[J]. Journal of Computer Applications, 2014, 34(12): 3605-3608.) |
[11] | 王益艳. 联合结构张量和变指数正则变分医学图像复原[J]. 计算机工程与应用, 2016, 52(15): 208-211. (WANG Y Y. Medical image restoration via joint structure tensor and variable index regularization variational[J]. Computer Engineering and Applications, 2016, 52(15): 208-211. DOI:10.3778/j.issn.1002-8331.1601-0331) |
[12] | KUIJPER A. p-Laplacian driven image processing[C]//Proceedings of the 2007 IEEE International Conference on Image Processing. Piscataway, NJ:IEEE, 2007, 5:V-257-V-260. |
[13] | HUNG K W, SIU W C. Robust soft-decision interpolation using weighted least squares[J]. IEEE Transactions on Image Processing, 2012, 21(3): 1061-1069. DOI:10.1109/TIP.2011.2168416 |
[14] | WANG Z, BOVIK A C, SHEIKH H R, et al. Image quality assessment:from error visibility to structural similarity[J]. IEEE Transactions on Image Processing, 2004, 13(4): 600-612. DOI:10.1109/TIP.2003.819861 |