摘要:本文提出一种偏微分方程约束的非局部均值图像插值模型。非局部均值模型在全局信息控制下处理图像，像素点灰度值为原始未插值图像全部像素的加权平均结果。偏微分方程在局部几何特征约束下处理图像，本文使用像素曲率偏微分方程作为非局部均值模型相似度函数的约束，使得与待处理像素具有相似几何特征的像素点具有较高的权重。线性特征位置的像素具有较高相似度，从而保持图像边缘特征使得处理后的图像清晰。传统模型插值放大后图像中存在严重的锯齿现象，本文使用各向异性扩散的偏微分方程对插值图像做预处理，沿边缘方向的前向扩散平滑锯齿。本文提出的新模型应用在预处理后的图像中平滑图像灰度块，保持图像边缘清晰，且最终结果无锯齿现象。理论和实验证明了本文提出的模型在图像插值中的有效性。
    关键词:图像插值　偏微分方程  非局部均值  锯齿现象

    1、引言
    在数字图像的传输过程中，为保证较高的传输速率，信道传输的图像分辨率低，而终端设备例如打印机显示器等分辨率高。因此插值放大提高图像分辨率是数字图像处理中一个非常重要的问题。传统插值模型例如最近邻，双线性和双三次模型通过匹配多项式函数得到放大后图像像素灰度值 [6]，像素重复显示使得图像中存在严重的锯齿现象。
    非局部均值 (Non-local means) 模型根据像素自相似性处理图像，插值放大后图像像素点灰度值为原始图像全部像素灰度值的加权平均 [7,8]。与待处理像素点相似的像素具有较高权重，相似度取决于像素邻域，而不是单个像素点灰度值。基于全局信息的非局部均值模型处理的图像平滑且保持图像中的纹理特征[7,8]。由于图像边缘像素具有较少相似像素点，因此非局部均值模型处理后图像边缘模糊 [1]。插值放大的图像中存在严重的锯齿现象，非局部均值模型不能消除锯齿：若模型中参数取较小值则图像仍存在锯齿现象不够平滑；若参数取较大值则图像过平滑，模糊边缘及纹理特征。
    偏微分方程(Partial differential equation, PDE)根据局部几何特征处理图像 [2-4]，基于各向异性扩散性能，PDE在平滑图像同时锐化边缘特征。PDE插值放大的图像清晰但被分割为邻接的灰度区域，因此图像不够平滑；PDE模型无法保持纹理特征，不能处理纹理图像。
    本文提出一种偏微分方程约束的非局部均值插值模型，新模型同时基于局部和全局信息处理图像。本文使用像素点曲率PDE约束邻域相似度函数，使得与待处理像素具有相似几何特征的像素具有较高权重，新模型保持图像线性特征，使得插值放大后的图像清晰。PDE具有各向异性扩散性能[2-5]，沿边缘方向前向扩散平滑锯齿，垂直边缘方向后向扩散锐化边缘。本文使用PDE预处理图像，去除插值图像中边缘位置的锯齿；PDE约束的非局部均值模型平滑预处理后的图像，保持图像中的边缘特征，插值放大的图像中不存在锯齿现象。
    本文由六部分组成，引言后，第二部分分析了偏微分方程(PDE）在图像插值中应用的原理，推导出插值预处理模型；第三部分分析本文提出的PDE约束的非局部均值插值模型；第四部分分析模型中PDE参数的数值实现格式；第五部分给出实验结果；最后给出一些结论，并提出一些有待继续研究的问题。
    2、PDE预处理
    图1显示PDE的两种扩散方向，沿水平曲线方向 和梯度方向 ， 垂直于 。沿水平曲线方向的扩散改变边缘位置像素灰度值

（1）

    其中x和y为像素欧式坐标， 为梯度模， 为水平曲线方向二阶导数， 为扩散系数。 前向扩散，像素灰度值增加相邻像素间的差异减小；  后向扩散，像素间差异放大。前向扩散平滑图像后向扩散锐化边缘。本文使用沿水平曲线方向的前向扩散平滑边缘的锯齿，即 。为了保持图像清晰，需要加入沿梯度方向的后向扩散锐化[3,4]。本文使用冲击滤波器方程[10]

                                                    图1 各向异性扩散方向
                                   Fig.1 Anisotropic diffusion directions （2）
    其中 为符号函数,  为沿梯度方向的二阶导数，G 为高斯滤波器，* 为卷积算子。高斯平滑后图像平滑区及锯齿边缘 为零，无后向扩散的锐化操作；边缘像素 不为0，后向扩散锐化边缘。本文使用的双向扩散PDE预处理模型为

(3）

    沿水平曲线方向前向扩散平滑图像中的锯齿现象，沿梯度方向的后向扩散锐化边缘。PDE预处理后的插值图像清晰且无锯齿效应。
    3、PDE约束的非局部均值插值模型
    图像中像素间关系符合马尔可夫性 [9]：任一像素灰度值仅与它邻域内像素有关，而与较远位置像素无关。两个像素具有相似的邻域特征则它们的几何特征和灰度值也相似。非局部均值模型对原始图像中全局信息加权平均得到待处理像素灰度值[1,7,8]。根据马尔可夫性权重应为两像素邻域间距离的函数，两个像素相似则它们的邻域距离较小。邻域距离计算公式为

（4）

    其中 为原始未插值图像， 为插值后图像, 表示邻域， 表示 中的邻域而 表示 中邻域,k 为插值放大倍数，i,j 为放大后图像中像素坐标l,m 为原始图像像素坐标。距离函数 具有多种形式：p=1表示Manhattan距离，p=2表示Euclidean距离。在图像插值中，邻域间距离用于衡量像素间相似度。若某像素邻域中像素与另一邻域中相应像素不相似，则该像素具有较低相似度，即邻域中的每一像素都影响相似度函数的结果。因此在非局部均值模型中，应当使用Euclidean距离，p=2。任意两像素间的差异通过平方运算放大。邻域距离公式为

(5)

    图像中曲线通过曲率驱动的法向向量运动生成[11]，曲率的PDE形式为

(6)

    其中u表示图像像素灰度值， 表示散度算子，图像中线性特征位置的像素点具有较大的曲率值。使用曲率约束相似度函数，可以使得与待修复像素具有相似曲率值即相似几何特征的像素点具有较高权重，这些像素点对邻域距离相似度函数影响较大，从而保持图像中的线性特征。下面给出曲率的相似度匹配函数

 (7)

    其中 为插值后图像局部曲率， 为原始未插值图像曲率，匹配函数仍使用Euclidean距离。曲率约束的邻域间距相似度函数为

(8)

   原始图像 中像素权重计算公式为

   (9)

 其中h为阈值。 中全部像素（l，m )灰度值与相应权重 相乘结果总和为非局部均值模型处理结果u(i.j) 的灰度值。为了简化描述，像素(i,j) 记作x，(l,m) 记作y 。非局部均值模型记作 

    (10)

    其中Z(X) 为归一化算子， 为原始图像域， 为插值后的图像域，u为最终结果。h 影响处理效果，通常取 ， 为 中像素的标准方差。h 值大，图像平滑；h 值小，图像清晰，当h 值取得很小处理后的图像几乎没有改变。
    在阈值h 取值相同的情况下，PDE约束的非局部均值模型根据局部模块的几何特征处理图像，保持图像中的线性特征，图像比无约束的模型处理结果清晰，但图像中仍存在锯齿现象。使用PDE对 做预处理，则 中无锯齿效应。在预处理后的图像 中使用新模型，图像清晰且无锯齿效应。由上述分析，PDE约束的非局部均值模型处理优点为：
    • 模型处理的图像平滑。平滑区域的像素值为原始图像中全局信息加权平均的结果；
    • 处理后的图像边缘清晰。PDE约束相似度函数，保持图像的边缘特征；
    • 图像无锯齿现象。边缘像素为PDE预处理的结果，沿水平曲线方向前向扩散平滑锯齿。
    4、 曲率PDE的数值实现
    PDE处理结果与所选用的差分格式有关。本文模型中相似度函数曲率约束项是散度算子的扩散结果，下面对它的差分格式进行分析。中心差分格式为二阶精度，而前向和后向差分格式为一阶精度。首先将散度算子作中心分：                   

其中，

，(11)

以图2中E方向为例。   在 点做中心差分：
 (12)

可以看到判断O点的扩散情况，用到了X点的值，超出了图像处理用到的像素9邻域范围 [12]。所以采用半点中心差分格式，即对e点做差分：

                                                   图2 半点差分像素结构
                                       Fig.2  Half-point differential structure       (13)
    这样对虚拟点的近似，使得图像处理为单个像素点的9邻域。

　　5、实验结果及分析
    本文使用Lena和Pepper局部图像验证PDE约束的非局部均值模型应用在图像插值放大中的有效性。最小均方误差估计(MSE)，峰值信噪比(PSNR)，信噪比(SNR)都为衡量插值效果的数学依据，通常认为这些数据都不能完全反应插值效果 [1,4]，有些方法以牺牲图像清晰度为代价得到较小MSE [3]，所以本文只对模型处理的图像进行主观比较。

　　图3为局部Lena图像。(a)为原始未插值图像。(b)为最近邻法插值处理的图像，帽沿部分存在严重的锯齿现象，脸部也不够平滑。(c)为非局部均值模型插值结果。Lena的脸部非常平滑，一些边缘特征例如帽沿和鼻子的轮廓不清晰，并且帽子和脸部边缘存在很严重的锯齿现象，非局部均值模型不能去除插值图像中的锯齿现象。(d)为本文提出的PDE约束的非局部均值模型处理结果。(d)比(c)清晰，图像不存在过平滑的现象，脸部的光照渐变特征都保留下来，脸部线性特征清晰，但图像中仍存在锯齿现象。(e)为节2中介绍的各向异性扩散的冲击滤波器PDE预处理后的图像。PDE根据几何特征局部处理图像，切向梯度方向的扩散改变图像对比度，(e)中图像不平滑，脸部存在分割的灰度块，但是图像中不存在锯齿现象，沿水平曲线方向的前向扩散平滑边缘位置的锯齿。(f)为PDE约束的非局部均值模型对PDE预处理后插值放大的结果。图像平滑且保留光照渐变特征，脸部线性特征清晰没有过平滑现象，图像中边缘特征连续，不存在锯齿现象。比较(b)-(f)的图像可以看出，PDE去除锯齿现象预处理后的图像应用本文提出的新模型插值可以得到最好的放大效果。

　　图4、5为局部Pepper图像的插值放大结果。处理结果与Lena图像一致，PDE约束的非局部均值模型保持线性特征，处理结果清晰。PDE预处理后的图像无锯齿现象，新模型应用在预处理后的图像中，可以得到平滑、清晰且无锯齿现象的插值放大图像。

    6、总 结
    本文提出一种PDE约束的非局部均值插值模型。非局部均值模型基于全局信息处理图像，插值放大的结果平滑，但模糊边缘特征。本文将像素曲率加入到相似度函数约束中，使得与待处理像素具有相同几何特征的像素具有较高权重，保持图像中边缘清晰。单纯采用非局部均值模型不能消除图像中的锯齿现象。本文采用各向异性扩散的PDE对插值图像做预处理，沿边缘方向的前向扩散平滑边缘位置的锯齿。新模型应用在预处理后的图像中保持边缘特征，平滑灰度块并且无锯齿效应，实验证明本文提出的新模型应用在预处理后的插值图像中具有很好的放大效果。
    非局部均值和PDE是两种广泛应用的图像插值模型，都存在一些问题。将两种模型结合同时基于局部和全局信息插值放大图像可以得到很好的效果。如何选取模型中的参数，如何选择具有完善数学理论依据的结合方式还需要进一步的研究，我们将在今后的文章中加以讨论。

   参考文献：

    [1] Wittman T. Mathematical techniques for image interpolation. Oral Exam Paper, July 2005
    [2] Caselles V, Morel J M, and Sbert C. An axiomatic approach to image interpolation[J]. IEEE.Trans. Image Processing, 1998, 7(3), 376~386.
    [3] Liu Z X, Wang H J, and Peng S L. An image magnification method using joint diffusion[J]. Journal of Computer Science and Technology, 2004, 19, 698~707.
    [4] Fu S J. Adaptive image interpolation using coupled bidirectional flow[A]. IEEE Int. Conf on Image Processing[C], 2005, 2, II~ 970~3.
    [5] Gilboa G, Sochen N, and Zeevi Y Y. Forward-and-backward diffusion processes for adaptive image enhancement and denoising[J]. IEEE. Trans. Image Processing, 2002,11(7), 689~703.
    [6] Hou H S, Andrews H C. Cubic splines for image interpolation and digital filtering[J]. IEEE Trans. Acoustics, Speech, Signal Processing, 1978, 26, 508~517.
    [7] Antoni B, Bartomeu C, Morel J M. On image denoising methods CMLA Preprint, CMLA, 2004
    [8] Antoni B, Bartomeu C, Morel J M. Denoising image sequences does not require motion estimation[A]. Proc. IEEE Conference on Advanced Video and Signal Based Surveillance[C], 2005, 70~74.
    [9] Zou Mouyan. Deconvolution and signal recovery[D], National Defence Industry Press, 2001. 邹谋炎。反卷积和信号复原[D].北京：国防工业出版社，2001.
    [10] Alvarez L, Mazorra L. Signal and image restoration using shock filters and anisotropic diffusion[A]. SIAM J. Numer. Anal.[C], 1994, 31(2),, 590~605.
    [11] Osher S, Sethian J A, Front propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations[J]. J Comput Phys, 1988, 79, 12~49. 
    [12] Chan T F, Osher S and Shen J. The digital TV filter and nonlinear denoising[J]. IEEE Trans. Image Processing, 2001, 10: 231~241