基于结构张量图像建模方法的滤波性能研究

http://www.paper.edu.cn - 1 - 基于结构张量图像建模方法的滤波性能研究1 邵文泽，韦志辉 南京理工大学计算机科学与技术学院，南京(210094) E-mail: shaowenze1010@yahoo.com.cn 摘 要: 本文关注基于结构张量的图像建模方法，对基于结构张量的偏微分方程(PDE)和变 分泛函方法的滤波性能进行了系统的分析研究。基于角形强度度量和水平线演化理论，首先 研究设计出一种用于角点增强的角形冲击滤波器，以克服边缘冲击滤波增强图像的不足。基 于边缘冲击滤波器和角形冲击滤波器，随后研究分析了基于扩散张量的各向异性PDE的滤波 性能。研究指出，Weickert提出的散度型各向异性PDE对应的是具有角点保持作用的平滑增 强滤波机制；而Tschumperlé提出的可计算迹型PDE是散度型各向异性PDE的一种退化情形， 对应的是平滑滤波机制且不具有角点保持作用。建立在上述工作基础上，研究分析了基于结 构张量变分泛函方法满足角点保持性的条件。最后，建立一种面向应用的统一正则PDE框架， 从而更直观、更有效地刻画在平坦区域、边缘结构、以及角形结构的滤波性能。图像滤波、 图像放大、超分辨率重建、图像修补等多种图像处理结果验证了统一正则PDE框架的有效性。 关键词：偏微分方程；冲击滤波；结构张量；角点保持；图像插值；图像修补；超分辨率 中图分类号：TP 391 1. 引言 自20世纪80年代末，数字图像处理的深入研究极大推动了基于Fourier分析方法的经典信 号处理学科的发展。经典的数字信号处理往往基于信号的平稳性或者高斯随机过程假设，因 此线性算法被认为是最优的。然而，实际的数字信号尤其是自然图像往往存在具有边缘、角 点等具有重要视觉意义的非平稳几何结构，很难用简单的高斯过程进行有效刻画，所以数字 图像处理逐渐倾向于非线性算法的研究。因此，数字图像处理的主流研究包括两个重要方面 [1]：一是数字图像的理解、表示与数学建模研究；二是基于图像模型的图像处理算法设计及 性能研究。有效的图像建模对后续图像处理算法的设计和性能具有至关重要的作用。经过近 20年的发展，研究者们已经从多个角度建立起图像的数学模型，主要包括[1,2]：基于马尔科 夫随机场(Markov random field; MRF)的统计图像建模方法，基于正则化 函数 excel方差函数excelsd函数已知函数     2 f x m x mx m      2 1 4 2拉格朗日函数pdf函数公式下载 空间的几何图像 建模方法，基于调和分析的几何多尺度图像建模方法，基于多重分形的几何图像建模方法， 以及基于人类视觉系统感知的图像建模方法。基于正则化函数空间的几何图像建模方法将导 致基于变分泛函和几何偏微分方程(partial differential equation; PDE)的图像处理方法，统称为 变分PDE方法。与其它方法相比，变分PDE方法在数学理论和数值计算上具有如下显著特点： 可以直接处理图像中视觉重要的几何特征，如梯度、切线、曲率等；可以有效地模拟具有视 觉意义的动态过程，如各向同性扩散、各向异性扩散等；可以运用成熟的数值计算方法或现 代智能优化算法，如模拟退火、遗传算法等。 本文重点关注基于结构张量(structure tensor)的图像建模方法[3-10]，对基于结构张量的偏 微分方程(PDE)和变分泛函方法的滤波性能进行了系统的分析研究。基于角形强度度量和水 平线演化理论，首先研究设计出一种用于角点增强的角形冲击滤波器，以克服边缘冲击滤波 1 Supported by the Natural Science Foundation of China under Grant No.60672074 (国家自然科学基金); by the National High-Tech Research and Development Plan of China under Grant No.2007AA12E100 (国家高技术研究 发展计划(863)); by the National Research Foundation for the Doctoral Program of Higher Education of China under Grant No.M200606018 (教育部博士点基金); by the Natural Science Foundation of Jiangsu Province of China under Grant No.BK2006569 (江苏省自然科学基金); and by the Science-Technology Creation Plan for Graduate Students of Jiangsu Province of China (江苏省高校研究生科技创新计划). http://www.paper.edu.cn - 2 - 增强图像的不足。基于边缘冲击滤波器和角形冲击滤波器，随后研究分析了基于扩散张量的 各向异性PDE的滤波性能。研究指出，Weickert提出的散度型各向异性PDE [3-7]对应的是具有 角点保持作用的平滑增强滤波机制；而Tschumperlé提出的可计算迹型PDE [8-10]是散度型各向 异性PDE的一种退化情形，对应的是平滑滤波机制且不具有角点保持作用。建立在上述 工作基础上，研究分析了基于结构张量变分泛函方法满足角点保持性的条件，并对文献中的 现有变分泛函进行了相应的验证讨论。最后，建立一种面向应用统一正则PDE框架，从而更 直观、更有效地刻画在平坦区域、边缘结构、以及角形结构的滤波性能。图像滤波、图像放 大、超分辨率重建、图像修补等多种图像处理结果验证了统一正则PDE框架的有效性。 本文第2部分对基于结构张量图像建模的相关方法及其性能分析作了简单回顾；第3部分 研究设计出一种用于角点增强的角形冲击滤波器；第4部分研究分析了基于扩散张量的各向 异性PDE的滤波性能，讨论了基于结构张量变分泛函方法具有角点保持性能的条件，并且建 立了一种面向应用的统一正则PDE框架；第5部分给出多种图像处理任务的实验结果；第6 部分给出结论。 2. 基于结构张量的图像建模方法回顾 2.1 基于结构张量的各向异性PDE 基于PDE的图像分析与处理思想最早出现在1983年Witkin发表的文章“尺度空间滤 波”[11]，指出信号与具有不同尺度(方差)高斯函数卷积等价于以信号为初值的热传导方程。 假设 0u 为初始信号，与 0u 相关的“尺度空间”分析归于求解下面的方程 ( ) ( ) ( )00u u x,t , u x, u x t ∂ = Δ =∂ (1) 然而，高斯函数卷积作为低通滤波对应的是线性运算，或者换句话说，热传导方程本质上对 应的是各向同性的扩散过程，图像中的边缘、角点等重要几何结构将被过分光滑，不利于后 续的图像处理操作，如边缘检测、图像分割。事实上，数字图像的非高斯过程要求图像处理 算法应该具有非线性和各向异性的特点，使得滤噪的同时能够保护图像中重要的几何信息。 基于PDE的非线性图像处理算法始于1987年著名的Perona-Malik各向异性扩散方程[12]的提 出，通过在热传导方程中引入取值在0和1之间的单调下降“边缘停止”函数 g ，实现边缘保持 的图像滤波 ( )( ) ( ) ( )00u g u u , u x, u xt∂ = ∇ ∇ ∇ =∂ (2) 虽然Perona-Malik方程能够一定程度上实现边缘保持的图像滤波，但仍存在严重缺陷，例如： 初始图像 0u 被噪声严重污染时，会导致大量伪边缘结构的出现；Perona-Malik方程本身是不 适定问 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 ，函数 g 的选取往往导致逆扩散造成解的不稳定性，事实上，Perona-Malik方程对 应一种平滑-增强型扩散机制[13]。尽管如此，Perona-Malik方程开创了基于PDE的非线性图像 处理的先河。之后的10多年时间里，基于PDE的多数图像处理研究都是基于Perona-Malik方 程展开的，包括1994年Weickert[5]最先提出的基于结构张量的散度型各向异性PDE和2002年 Tschumperlé[8,9]最先提出的面向应用的迹型各向异性PDE。 首先引入结构张量的概念[5]。结构张量定义为： , 1,2; 1,2( ) ( ) ( )m n m nu G u u Jρ σ ρ σ σ = =∇ = ∗ ∇ ⊗∇ =J ， 其中 uσ 为图像 u 经过高斯滤波的图像(方差 0σ > )，减少噪声对求导运算的影响，从而增强边 缘检测的噪声鲁棒性；对张量 u uσ σ∇ ⊗∇ 进行高斯滤波(方差 0ρ ≥ )，以增强几何结构方向的 http://www.paper.edu.cn - 3 - 估计鲁棒性。矩阵 ρJ 对称且半正定，存在正交单位特征向量，分别记为 w 和 ⊥= wν 。其中， ( )T2 212 22 11 22 11 122 , ( ) 4 , / | |J J J J J J= − + − + =w w w w 指向几何结构的最大对比度方向，称为主向量；相应地，ν 指向几何结构的最小对比度方向， 看 作 是 几 何 结 构 的 方 向 。 它 们 相 应 的 特 征 值 分 别 记 为 μ 和 μ⊥ ， 定 义 为 2 2 0.5 11 22 22 11 12( (( ) 4 ) )J J J J Jμ = + + − + ， 2 2 0.511 22 22 11 12( (( ) 4 ) )J J J J Jμ⊥ = + − − + 。这两个值可以作为局部 几何结构的属性描述子：在平滑区域， 0μ μ⊥≈ ≈ ；在边缘区域， 0μ μ⊥>> ≈ ；在角形区域， 0μ μ⊥≥ >> 。 Weickert基于结构张量提出如下散度型各向异性PDE ( ) ( ) ( )0( ( )) 0u D u u , u x, u x t ρ σ ∂ = ∇ ∇ ⋅∇ =∂ J (3) 其中， ( ( ))D uρ σ∇J 称为扩散张量，对应2×2的(半)正定矩阵。显然，式(3)是Perona-Malik方程 (2)的推广，当 D=g ⋅ I 时式(3)即退化为式(2)，其中I是单位矩阵。一方面，Weickert从理论上 讨论了式(3)解的存在性、唯一性、极值性等问题；另一方面，Weickert等人将式(3)广泛应用 于图像去噪、指纹增强、光流场计算等多种图像处理问题中[3,4,14,15]。在图像处理的实际应用 中，Weickert将扩散张量定义为 T T1 2( ( )) =D uρ σ λ λ∇ +J ww νν ，根据特征向量 ,νw 确定PDE的主 要扩散方向，根据 1 2,λ λ 确定两个扩散方向上的扩散强度。因此，PDE式(3)的滤波性能主要 由扩散张量 ( ( ))D uρ σ∇J 具体刻画。 王正明等人[16]在2006年出版的《SAR图像提高分辨率技术》一书中关于PDE式(3)的滤 波性能给出了如下定理。 定理[16] 设 0u 为原始图像， T T1 2( ( )) =D uρ σ λ λ∇ +J ww νν 为连续的2×2扩散矩阵， ,νw 为其 规范 编程规范下载gsp规范下载钢格栅规范下载警徽规范下载建设厅规范下载 正交的特征向量分别对应图像的梯度方向和边缘方向， 1 2,λ λ 为其相应的特征根，利用扩散 平滑模型(3)对 0u 进行光滑近似于分别以 1 2,λ λ 的速度沿 ,νw 方向光滑。 Albert等人[13]在2002年出版的《Mathematical Problems in Image Processing》一书讨论基 于PDE的图像处理方法部分中，将主要PDE分成三类：平滑型(smoothing)PDE，平滑-增强型 (smoothing-enhancing) PDE，增强型(enhancing) PDE，并将PDE式(3)划分为平滑型PDE。 Tschumperlé等人[8-10]在2002年-2005年先后发表了多篇学术 论文 政研论文下载论文大学下载论文大学下载关于长拳的论文浙大论文封面下载 研究讨论了PDE式(3)的 滤波性能，指出扩散张量驱动的散度型PDE不能够完全刻画扩散张量所要体现的滤波性能， 并且提出一种扩散张量驱动的可计算迹型PDE ( ) ( ) ( )0t race ( ( )) 0u D u H , u x, u x t ρ σ ∂ = ∇ ⋅ =∂ J (4) 其中， H 为Hessian矩阵。Tschumperlé等人将PDE式(4)广泛应用于图像去噪、图像修补、图 像放大等多种图像处理问题中。 事实上，本文研究分析将指出，PDE式(3)实质上对应的是平滑-增强型(smoothing- enhancing)滤波机制，而PDE式(4)对应的才是平滑型(smoothing)滤波机制。利用PDE式(3)对 图像处理时，PDE不仅分别以 1 2,λ λ 的速度沿 ,νw 方向进行平滑，同时还具有增强角点和强边 缘结构的作用。 2.2 基于结构张量的变分泛函 基于正则化函数空间的几何图像建模方法将导致基于变分泛函和几何PDE的图像处理 方法。与2.1节中的边缘保持PDE类似，研究者们相应地提出如下具有边缘保持性能的变分 http://www.paper.edu.cn - 4 - 泛函[17] ( ) ( )rE u u dxϕΩ= ∇∫ (5) 其中， ( )ϕ ⋅ 称为势函数，满足如下基本条件： ( ) 0xϕ ≥ ， ( ) = (- )x xϕ ϕ ， ( ) 0' xϕ ≥ ，并且 ( )/' x xϕ 在 (0, )∞ 上是单调下降函数。例如，当 ( ) =| |x xϕ 时，式(3)即为TV(total variation)模型[18,19]；当 ( )xϕ 取为Huber范数时，式(3)即为Huber-MRF模型[20]。对应于变分泛函(5)的PDE为 ( )' uu u t u ϕ⎛ ⎞∇∂ = ∇ ⋅ ⋅∇⎜ ⎟⎜ ⎟∂ ∇⎝ ⎠ (6) 2002年，Weickert[21]首次提出如下基于结构张量的变分泛函用于运动估计 ( ) ( ) ( )( )T TtracerE u dϕ μ ϕ μ⊥Ω= + Ω∫ ww νν (7) 其中， ( )ϕ ⋅ 取为上述势函数。对应于式(7)的PDE为 ( ) ( )( )( )T Tu ' ' ut ϕ μ ϕ μ ⊥∂ = ∇ + ⋅∇∂ ww νν (8) 2005年，Tschumperlé等人[10]基于结构张量提出更为一般的变分泛函 ( ) ( ),rE u dψ μ μ⊥Ω= Ω∫ (9) 其中， ( ),ψ ⋅ ⋅ 通常采用势函数 ( )ϕ ⋅ 定义。当结构张量 ( )uρ σ∇J 中的 ρ 和σ 均取0时，式(9)即退化 为式(5)。对应于式(9)的PDE为 2 T Tu u t ψ ψ μ μ ⊥ ⎛ ⎞⎛ ⎞∂ ∂ ∂= ⋅∇ + ⋅∇⎜ ⎟⎜ ⎟∂ ∂ ∂⎝ ⎠⎝ ⎠ ννww (10) 事实上，比较PDE式(8)和PDE式(9)知，式(7)也是式(9)的一种特殊情形。 注意到，Tschumperlé等人在提出变分泛函式(9)时没有给出 ( ),ψ ⋅ ⋅ 的选取条件。事实上， 不同的 ( ),ψ ⋅ ⋅ 选取方式对于式(9)滤波性能将有不同的影响，本文重点研究分析了结构张量驱 动的变分泛函式(9)满足角点保持性的条件，并对文献中的现有变分泛函进行了验证讨论。 3. 角形冲击滤波器 3.1 边缘冲击滤波器 引入角形冲击滤波之前，首先介绍边缘冲击滤波器[22]。 在图像处理中，增强型(enhancing)PDE主要用于增强图像从而提高图像的清晰度，对应 着图像平滑(smoothing)的逆过程。Osher等人 [22]在1992年首次提出一种称为冲击滤波器 (shock filter)的PDE图像增强方法。最常用的冲击滤波器定义如下 ( )( )2u sign D , u t u η η∂ = − ∇∂ (11) 其中， / || ||u u= ∇ ∇η 为图像的梯度方向。可见，冲击滤波器根据边缘检测项 2 ( )D ,u η η 的符号 确定冲击流的方向，根据图像的梯度模确定冲击流的强度。然而，PDE(11)存在两点重要不 足：其一，增强边缘的同时增强噪声点；其二，只能增强边缘不能增强角形结构，因此实 质上是边缘冲击滤波。针对第一点不足，研究者们通过耦合边缘冲击滤波和平滑型PDE[23]， 实现同时增强边缘和抑制噪声的目的。然而，对于如何在PDE框架下实现角形增强，至今没 有解决方案。为此，对应于边缘冲击滤波，本文研究设计出一种角形冲击滤波器，通过耦 合边缘冲击滤波和角形冲击滤波，能够同时实现增强边缘和角点结构的目的。 受边缘冲击滤波器的启发，角形冲击滤波器的设计分成两步：第一，通过研究角形强 度度量确定角形冲击流的速度；第二，在基于水平集方法的曲线演化理论框架下研究如何 http://www.paper.edu.cn - 5 - 确定角形冲击流的方向。 3.2 角形强度度量(MoCS; measures of corner strength) 由于角形强度是用于确定角形冲击流的速度，因此本文认为角形强度度量应该满足以下 三个基本 准则 租赁准则应用指南下载租赁准则应用指南下载租赁准则应用指南下载租赁准则应用指南下载租赁准则应用指南下载 ：噪声鲁棒性；边缘鲁棒性；度量精确性。上述准则的提出主要基于以下考虑： (1) 角形的强度度量应该对噪声鲁棒，角形冲击滤波器不能同时增强噪声； (2) 角形的强度度量只能刻画角形区域的强度，不能刻画其它结构的强度； (3) 传统的角点检测算法不能替代角形的强度度量，只能检测锐化的角点； (4) 不同的角形具有不同的强度，对应的角形冲击流应该具有不同的速度。 目前，图像处理文献中对于角形强度度量的研究还相对较少，主要有以下两种度量[24]： 2 2 1 2 2 2xx y xy x y yy x x y MoCS u u u u u u u u u u u κ− += = = ⋅ ∇+ ξξ (12) ( ) ( )2 22 1 1 ( ) /( )MoCS moa u uμ μ μ μ⊥ ⊥= − ⋅ ∇ = − − + ⋅ ∇ (13) 其中， ⊥ξ η ， κ 为曲率， 2 2( ) /( )moa μ μ μ μ⊥ ⊥= − + 为各向异性测度。由于结构张量的特征值 可以作为局部几何结构的属性描述子，即在平滑区域， 0μ μ⊥≈ ≈ ；在边缘区域， 0μ μ⊥>> ≈ ； 在角形区域， 0μ μ⊥≥ >> ，因此 moa 实际上起到指示边缘区域的作用。然而，式(12)和式(13) 对噪声过分敏感，另外式(12)对曲线型边缘较为敏感，因此需要研究设计更为精确鲁棒的角 形强度度量。 首先，定义 2×2 矩阵的散度运算。对于矩阵 1 2; 1 2( )i , j i , j ,M m = == ，定义 T 11 12 21 22m m m m M , x y x y ⎛ ⎞∂ ∂ ∂ ∂∇ ⋅ = + +⎜ ⎟∂ ∂ ∂ ∂⎝ ⎠ (14) 在引入新的角形强度度量之前，我们给出如下定理。 定 理 1 假 设 ( , )u x y 处 处 二 阶 连 续 可 微 ， ξ 为 图 像 的 切 线 方 向 ， 则 等 式 ( )( ) ( )TT 2 ,u D u∇ ⋅ ⋅∇ = −ξξ ξ ξ 成立。 证明：由 ( , )u x y 的二阶连续可微性知， xy yxu u= 成立。根据矩阵的散度运算定义，有 ( ) ( )( ) ( ) ( )( ) ( ) 2 2 2 T 2 2 2 2 22 4 2 2 2 + 21 2 y x y x y x y yx y x xx y yx xy y x yy x y x xy y yy xx y x xy x y x xx u u u x yu u u u u x yu u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u ⎛ ⎞⎛ ⎞ ⎛ ⎞−∂ ∂⎜ ⎟⎜ ⎟ ⎜ ⎟+⎜ ⎟ ⎜ ⎟∂ ∂⎜ ⎟∇ ∇⎝ ⎠ ⎝ ⎠⎜ ⎟∇ ⋅ = ⎜ ⎟⎛ ⎞ ⎛ ⎞−∂ ∂⎜ ⎟⎜ ⎟ ⎜ ⎟+⎜ ⎟ ⎜ ⎟⎜ ⎟∂ ∂∇ ∇⎝ ⎠ ⎝ ⎠⎝ ⎠ ∇ − + + − − ∇ + = ∇ − − ∇ + + ξξ ( )( ) ( )( )2 2 2 2 2 3 3 4 3 3 2 2 2 2 - 2 3 21 3 2 y xy x xy x x xy y yy x y xy x y yy x y xx xy y x yy x xy xx y x y xx x y xy x y yy u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u uu ⎛ ⎞⎜ ⎟⎜ ⎟+ ∇ +⎜ ⎟⎝ ⎠ + − − −⎛ ⎞= ⎜ ⎟− − + + −∇ ⎝ ⎠ 从而 http://www.paper.edu.cn - 6 - ( )( ) ( ) ( ) ( )( ) T T 3 2 2 2 2 3 4 4 4 2 2 2 2 2 2 4 2 2 1 2 2 1 + 2 + 2 + 2 x y xy x y yy x y xx x y xy x yy y xx x y xx x y xy x yy y y xx x y xy x yy y xx x y xy x y u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u u ∇ ⋅ ⋅∇ = − − + − −∇ = − − + − −∇ − −= ξξ ( ) 2 2 2 y x yu u D u , + = − ξ ξ ■ 根据结构张量的定义，当 ( )uρ σ∇J 中的 ρ 和σ 均取 0 时， ( )uρ σ∇J 的特征向量 ,νw 便分别 退化为 ,η ξ 。由于结构张量通过考虑局部邻域增强结构方向估计的鲁棒性，更重要的是，只 有考虑局部邻域才能识别角形结构，因此 ,w ν 比起 ,η ξ 能够更为精确地刻画局部结构的方向。 为此，根据 1MoCS 和定理 1，我们提出如下新型角形强度度量 ( )( )TT3MoCS u= ∇ ⋅ ⋅∇νν (15) 图 1 和图 2 给出两组实验比较三种角形强度度量检测角形区域的有效性，分别对应无噪情形 和噪声情形。在结构张量中， ρ 和 σ 分别取值 2 和 1.5。角形强度度量图像利用公式 5 50, = 1 , 2 ,3jMoCS j ∗ + 进行计算显示。对于两种情形下的角形度量图像，通过综合考虑噪声 鲁棒性、边缘鲁棒性、以及度量精确性，本文认为 3MoCS 比 1MoCS 和 2MoCS 更适于作为角 形冲击滤波的冲击速度。 图 1 不含噪声情形. 从左至右分别为原始图像和利用 MoCS1、MoCS2、MoCS3得到的角形强度度量图像 Fig. 1 Noiseless case. From left to right: original image and measured images calculated respectively by MoCS1、 MoCS2 and MoCS3 图 2 高斯噪声情形. 从左至右分别为噪声图像和利用 MoCS1、MoCS2、MoCS3得到的角形强度度量图像 Fig. 2 Gaussian noise case. From left to right: original image and measured images calculated respectively by MoCS1、MoCS2 and MoCS3 3.3 角形冲击滤波器 角形冲击滤波器设计的第二步是研究如何确定角形冲击流的方向。事实上，本文角形冲 击流的方向是在基于水平集方法的曲线演化理论框架下确定的。首先，我们给出角形冲击滤 波器的一般形式 http://www.paper.edu.cn - 7 - ( )( ) ( )2 ,u sign D MoCS t u uξ ξ∂ = − ⋅∂ (16) 其中， ( ) 0MoCS u ≥ 为某种角形强度度量，具备噪声鲁棒性、边缘鲁棒性、以及度量精确性。 为了说明 PDE 式(16)具有增强角形区域的作用，我们暂且考虑如下 PDE ( )( ) ( )2 ,u sign D MoCS t u uξ ξ∂ = ⋅∂ (17) 注意到 ( ) 0MoCS u ≥ 作为角形强度度量具备噪声鲁棒性、边缘鲁棒性、以及度量精确性，式 (17)实际上只作用于角形区域。由于图像可以利用其水平集来表示，因此式(17)可以基于曲 线演化理论进行分析。根据水平线演化理论[13,18,25]，给定初始曲线 0 ( )pC ，可以通过定义一 个高维函数 2( , , ) : [0, )x y t R T Rφ × → 来表示： 0 2( ) = {( ) : ( 0) = 0}, , ,p Rx y x yφ∈C 。例如， ( 0), ,x yφ 可 定义为 ( 0) =, ,x y dφ ± ，其中 d 为点(x, y)到曲线 0 ( )pC 的距离。考虑如下演化 PDE ( )( ))sign MoCS t φ κ φ φ∂ = ∇ ⋅∂ (18) 其中， ( 0), ,x yφ 为初始条件。根据曲率驱动的水平线演化理论知[13,18,25]，PDE 式(18)演化的 速度大小为 ( )/ || ||cm φ φ∇ ，演化的方向由曲率 κ 的符号确定，并且演化曲线始终对应函数 ( ), ,x y tφ 的零水平集。当 d 在曲线内部取负而在曲线外部取正，对于凸角形区域κ 取正，此 时角形向内部演化；而对于凹角形区域κ 取负，此时角形向外部演化。相反地，当 d 在曲线 外部取负而在曲线内部取正，对于凸角形区域κ 取负，此时角形向外部演化；而对于凹角形 区域κ 取正，此时角形向内部演化。因此，PDE 式(17)实际上起到平滑角形结构的作用，从 而得到圆滑的角形结构。由于 PDE 式(17)是 PDE 式(16)的反过程，所以 PDE 式(16)确实具 有增强角形区域的作用。当角形强度度量取为 3MoCS 时，角形冲击滤波器即为 ( )( ) ( )( )T2 T,u sign D u t u ξ ξ∂ = − ⋅ ∇ ⋅ ⋅∇∂ νν (19) 图 3 分别给出无噪情形和含噪情形下 PDE 式(19)增强角形区域的效果图。可见，本文角 形冲击滤波器不仅能够较好地增强角形结构，同时具有较强的噪声鲁棒性。图 4 给出了高斯 模糊情形下基于边缘冲击滤波器和角形冲击滤波的增强效果，较好地显示了本文角形冲击滤 波器研究的必要性和重要性。 图 3 角形冲击滤波. 从左至右分别为原始图像及相应的角形增强图像、噪声图像以及相应的角形增强图像 Fig. 3 Corner shock filtering. From left to right: original image and the corresponding corner-enhanced image, the noisy image and the corresponding corner-enhanced image http://www.paper.edu.cn - 8 - 图 4 图像增强. 从左至右分别为高斯模糊图像、基于边缘冲击滤波的增强图像、基于角形冲击滤波的 增强图像、以及通过耦合边缘冲击滤波和角形冲击滤波得到的增强图像 Fig. 4 Image enhancement. From left to right: Gaussian-blurred image, enhanced image by edge shock filtering, enhanced image by corner shock filtering, and enhanced image by coupling the edge shock filtering and the corner shock filtering 4. 基于结构张量的图像建模滤波性能分析及统一正则 PDE 框架 4.1 基于结构张量的各向异性PDE滤波性能分析 建立在上述工作基础上，本节具体研究分析基于结构张量的各向异性PDE的滤波性能。 我们首先给出如下定理。 定理2 对于扩散张量 T T0 0 1 2( ( )) = (| | ) ( | | )D u u | | u | |λ λ∇ ∇ + ∇J ηη ξξ ，规范正交的特征向量 ,η ξ 分别 对应图像的梯度方向和切线方向， 1 2,λ λ 分别对应特征向量的两个特征根，则下述等价式成 立 ( )T T1 2 1( (|| ||) (|| |) ) ( (|| ||) )u uu u| u u ut tλ λ λ∂ ∂= ∇ ⋅ ∇ + ∇ ⋅∇ ⇔ = ∇ ⋅ ∇ ⋅ ∇∂ ∂ηη ξξ 证明：事实上， 1T2 T (( (|| |) ) ) + (( (|| ||) )u / t u | u u uλ λ∂ ∂ = ∇ ⋅ ∇ ⋅∇ ∇ ⋅ ∇ ⋅ ∇ξξ ηη 2 T 2 2 T 1 1 2 2 2 2 (|| ||) ( ) ( ) [ ( (|| ||) )] + (|| ||) ( ) ( ) [ ( (|| ||) )] = (|| ||) ( ) - (|| | u / t u D u u u u D u u u u D u u λ λ λ λ λ λ ⇔ ∂ ∂ = ∇ ⋅ ⊗ + ∇ ⋅ ∇ ⋅ ∇ ⋅ ⊗ ∇ ⋅ ⊗ + ∇ ⋅ ∇ ⋅ ∇ ⋅ ⊗ ∇ ⋅ ⊗ ∇ ξ ξ ξ ξ η η η η ξ ξ 2 2 2 T 1 1 1 T 1 1 1 |) ( ) + (|| |) ( ) + (|| ||) ( ) + ( ) (|| ||) = (|| ||) ( ) (|| ||) = (|| | D u u | D u u D u u u u u + u u u λ λ λ λ λ λ ⋅ ⊗ ∇ ⋅ ⊗ ∇ ⋅ ⊗ ∇ ⋅∇ ∇ ∇ ⋅ Δ ∇ ⋅∇ ∇ ∇ ξ ξ η η ξ ξ 2 2 1 1 2 2 1 1 1 |) (|| ||) ( ) + [|| || (|| ||)]' ( ) = (|| ||) ( ) + [|| || (|| ||)]' ( ) = ( (|| ||) ) u - u D u u u D u u D u u u D u u u λ λ λ λ λ ⋅ Δ ∇ ⋅ ⊗ ∇ ⋅ ∇ ⋅ ⊗ ∇ ⋅ ⊗ ∇ ⋅ ∇ ⋅ ⊗ ∇ ⋅ ∇ ⋅∇ η η η η ξ ξ η η ■ 由定理2知，扩散张量 T T0 0 1 2( ( )) = (| | ) ( | | )D u u | | u | |λ λ∇ ∇ + ∇J ηη ξξ 驱动的散度型各向异性PDE式 (3)实质上等价于各向异性Perona-Malik方程，扩散强度仅由η 方向上的扩散系数 1( | | )u | |λ ∇ 决 定。可见，基于扩散张量的散度型PDE不能较好地刻画扩散张量所要表达的滤波性能。因此， 王正明等人[16]关于PDE式(3)的滤波性能所给出的定理是不够准确的。此外注意到，此时PDE 式(3)是一种平滑-增强型(smoothing-enhancing)PDE，具有平滑-增强的滤波性能，Albert等人 [13]将PDE式(3)划分为平滑型(smoothing) PDE也是不合适的。 对于Tschumperlé等人[8-10]提出的基于扩散张量的可计算迹型PDE式(4)，我们给出如下定 理。 http://www.paper.edu.cn - 9 - 定理3 假设 ( , )u x y 处处二阶连续可微， T T1 2( ( )) =D uρ σ λ λ∇ +J ww νν 为连续的2×2扩散矩阵，规 范正交的特征向量 ,ν w 分别对应图像的边缘方向和垂直方向， 2 1,λ λ 分别对应特征向量的两个 特征根，利用可计算迹型PDE式(4)扩散近似于分别以 1 2,λ λ 的速度沿 ,νw 方向光滑。 证明：根据迹及Hessian矩阵的定义，有 ( )T T1 2 T T 1 2 T T1 2 2 11 1211 1 211 12 22 t race ( ) t race ( ) t race trace xx xy yx yy xx xy xx xy yx yy yx yy u / t H u u u / t u u u u u u u u u u λ λ λ λ λ λ λ ∂ ∂ = + ⋅ ⎛ ⎞⎛ ⎞⇔ ∂ ∂ = + ⋅⎜ ⎟⎜ ⎟⎝ ⎠⎝ ⎠ ⎛ ⎞⎛ ⎞ ⎛ ⎞= ⋅ + ⋅⎜ ⎟⎜ ⎟ ⎜ ⎟⎝ ⎠ ⎝ ⎠⎝ ⎠ ⎛ ⎞= ⋅⎜ ⎟⎝ ⎠ νν νν νν ww ww ww w w w w w w ( ) ( ) 2 2 2 11 1211 2 211 12 22 2 2 2 21 11 12 2 11 1211 22 11 22 1 2 + 2 + + 2 + ( ) ( ) xx xy xx xy yx yy yx yy xx xy yy xx xy yy u u u u u u u u u u u u u u D u D u λ λ λ λ λ ⎛ ⎞⎛ ⎞ ⎛ ⎞⎛ ⎞+ ⋅⎜ ⎟⎜ ⎟ ⎜ ⎟⎜ ⎟⎝ ⎠⎝ ⎠ ⎝ ⎠⎝ ⎠ = + = ⊗ + ⊗ ν ν ν ν ν ν ν ν ν ν ν ν w w w w w w 其中， T T11 12 11 12( , ) ( , ),= =ν ν νw w w 。 ■ 可见，基于扩散张量 T T1 2( ( )) =D uρ σ λ λ∇ +J ww νν 的迹型各向异性PDE式(4)能够较好地刻画扩 散张量所要表达的滤波性能。同时我们注意到，式(4)是一种平滑型(smoothing) PDE。然而， 下文研究将指出：可计算迹型PDE式(4)是散度型各向异性PDE式(3)的一种退化情形，式(4) 虽然能够直观地刻画PDE在图像平坦区域以及边缘区域的滤波性能，但是没有刻画角形区域 的滤波性能，使得实际图像处理过程中角点结构往往被破坏。 考虑一般情形下扩散张量 T T1 2( ( )) =D uρ σ λ λ∇ +J ww νν ( 0, 0σ ρ> > )驱动的PDE式(3)的滤 波性能。之前，我们首先给出如下定理。 定理4 对于规范正交向量组 w 和ν ，等式 ( ) ( )( ) ( ) ( )( )T TT Tu u∇ ⋅ ∇ ⋅ = − ∇ ⋅ ∇ ⋅ ννww 成立。 证明：令 (cos , s in )θ θ=ν ，则 ( s in , cos )θ θ= −w 。根据矩阵散度运算法则，得 ( )( ) ( ) T 2T T 2 T T 2 2 2 cos cos sin sin cos sin cos cos sin 1 0 sin cos sin 0 1 cos 1 cos u u u u θ θ θ θ θ θ θ θ θ θ θ θ θ θ ⎛ ⎞⎛ ⎞− ∇ ⋅ ⋅∇ = − ∇ ⋅ ⋅∇⎜ ⎟⎜ ⎟⎝ ⎠⎝ ⎠ ⎛ ⎞ ⎛ ⎞− −⎛ ⎞ ⎛ ⎞= ∇ ⋅ ⋅∇ + ∇ ⋅ ⋅∇⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟− −⎝ ⎠ ⎝ ⎠⎝ ⎠ ⎝ ⎠ − − −= ∇ ⋅ νν ( ) ( )( ) T 2 T 2 T T 2 sin sin cos sin 1 sin cos sin sin cos cos u u u θ θ θ θ θ θ θ θ θ θ ⎛ ⎞⎛ ⎞ ⋅∇⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟− − −⎝ ⎠⎝ ⎠ ⎛ ⎞−⎛ ⎞= ∇ ⋅ ⋅∇ = ∇ ⋅ ⋅∇⎜ ⎟⎜ ⎟−⎝ ⎠⎝ ⎠ ww ■ 关于扩散张量 T T1 2( ( )) =D uρ σ λ λ∇ +J ww νν ( 0, 0σ ρ> > )驱动的PDE式(3)的滤波性能，我们 以定理的形式给出如下描述。大量的文献查阅表明，本文定理的证明过程及结论在国内及国 际尚属首创。 定理5 假设 ( , )u x y 处处二阶连续可微， T T1 2( ( )) =D uρ σ λ λ∇ +J ww νν ( 0, 0σ ρ> > )为连续的2×2扩 散矩阵，规范正交的特征向量 ,ν w 分别对应图像的边缘方向和垂直方向， 2 1,λ λ 分别对应特征 向量的两个特征根，利用各向异性PDE式(3)进行扩散时，PDE不仅具有以 1 2,λ λ 为速度沿 ,νw 方向光滑的作用，同时还具有增强角点和强边缘结构的作用。因此，各向异性PDE式(3)对应 的是平滑-增强型(smoothing-enhancing)滤波机制。 证明：首先对PDE式(3)作如下分解： http://www.paper.edu.cn - 10 - 由分解结果知，第 1 项正是Tschumperlé等人提出的可计算迹型PDE式(4)。因此，利用各向 异性PDE式(3)进行扩散时，PDE具有以 1 2,λ λ 为速度沿 ,νw 方向光滑的作用。 下面，具体分析第 2 项和第 3 项在滤波过程中的作用。根据定理4，有 首 先分析第 4 项。为了实现边缘保持性或边缘增强性的图像处理任务(如图像去噪、图像放大 等)，在图像边缘、角点等几何结构位置， 1 2,λ λ 之间应该满足如下关系： 2 1λ λ> 。另一方面， 注意到规范正交向量 ,ν w 是结构张量 ( )uρ σ∇J ( 0, 0σ ρ> > )的特征向量，规范正交向量 ,η ξ 是 结构张量 0 0( )u∇J 的特征向量，而 ,ν w 和 ,η ξ 两者都是用来刻画图像局部几何结构方向的，只 是 ,ν w 比起 ,η ξ 相对而言具有较高的精度。我们不难得出结论：在 0, 0σ ρ> > 取值相对较小时， 下式成立 ( )( )( ) ( )( )( ) ( )( )T TT T 2sign u sign u sign D u ,∇ ⋅ ⋅∇ = ∇ ⋅ ⋅∇ = −νν ξξ ξ ξ 从而得到 ( ) ( ) ( )( ) ( ) ( )( )( ) ( ) ( )( ) ( ) ( )( )( ) ( ) ( )( ) ( )( ) ( ) ( )( ) T T TT T T2 1 T TT T T2 T 4 u s ign u u s ign u u sign D u , u λ λ− = ∇ ⋅ ∇ ⋅ = ∇ ⋅ ∇ ⋅ ⋅ ∇ ⋅ ∇ ⋅ = ∇ ⋅ ∇ ⋅ ⋅ ∇ ⋅ ∇ ⋅ = − ⋅ ∇ ⋅ ∇ ⋅ νν νν νν ξξ νν ξ ξ νν 在角点等几何结构位置， 1 2,λ λ 之间应该满足关系 2 1λ λ> ，因此根据第3部分关于角形冲击滤 波器的研究，我们得出结论：第 4 项起到增强角点增强的作用。 其次分析第 5 项。为了实现边缘保持或边缘增强性的图像处理任务(如图像去噪、图像 放大等)， 1 2,λ λ 之间除了满足关系 2 1λ λ> 外， 1λ 通常选取为单调下降函数， 2λ 通常选取为单 调下降函数或常值函数。特别地，当 2λ 取为常值函数时， ( ) ( ) ( )T T 14 u λ= ∇ ⋅ ⎡ ⋅ ∇ ⎤⎣ ⎦ww 因此，我们考虑 1 2,λ λ 均为单调下降函数的情形。根据关于Perona-Malik方程的已知研究，当 1 2,λ λ 均为单调下降函数时， ( )T T1 2 2 211 12 11 1211 11 1 22 211 12 11 1222 22 2 21 1 11 12 2 2 11 1211 11 2 21 1 11 12 2 2 11 1222 22 ( ) ( ) x y x y y x y x u / t u u / t u u u u u u u u u λ λ λ λ λ λ λ λ λ λ λ λ ∂ ∂ = ∇ + ⋅∇ ⇔ ⎛ ⎞⎛ ⎞ ⎛ ⎞∂ ∂ = ∇ + ⋅∇⎜ ⎟⎜ ⎟ ⎜ ⎟⎝ ⎠ ⎝ ⎠⎝ ⎠ ⎛ + +⎛ ⎞ ⎛ ⎞= ∇ +⎜ ⎟ ⎜ ⎟+ +⎝ ⎠ ⎝ ⎠⎝ νν ν ν ν ν ν ν ν ν ν ν ν ν ww w w w w w w w w w w w w ( ) ( ) ( ) ( )( ) ( ) ( )( ) ( ) ( )( ) ( ) ( )( ) ( ) ( ) 2 2 T T2 2 2 2 T T1 11 12 2 11 12 1 211 22 11 22 T TT T1 2 1 2 TT T 1 2 + 2 + + 2 + ( ) ( ) + trace ( ) + xx xy yy xx xy yy u u u u u u u u D u D u u u H u λ λ λ λ λ λ λ λ λ λ λ ⎞⎜ ⎟⎠ ⎡ ⎤ ⎡ ⎤= + + ∇ ⋅∇ ⋅ + ∇ ⋅∇ ⋅⎣ ⎦ ⎣ ⎦ ⎡ ⎤ ⎡ ⎤= ⊗ + ⊗ ∇ ⋅∇ ⋅ + ∇ ⋅∇ ⋅⎣ ⎦ ⎣ ⎦ = + ⋅ ∇ ⋅∇ ν ν ν ν νν ν ν νν νν w w w w ww w w ww ww ( )( ) ( ) ( )( )TT T1 2 1 2 3 u λ⎡ ⎤ ⎡ ⎤⋅ + ∇ ⋅∇ ⋅⎣ ⎦ ⎣ ⎦ + + � ννww ( ) ( )( ) ( ) ( )( ) ( ) ( )( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) T TT T1 2 T T T TT T T T 1 2 1 2 T TT T T 2 1 1 2 2 3 4 5 u u u u u u u u λ λ λ λ λ λ λ λ λ λ Δ + ⎡ ⎤ ⎡ ⎤= ∇ ⋅∇ ⋅ + ∇ ⋅∇ ⋅⎣ ⎦ ⎣ ⎦ = ⋅ ∇ ⋅ ∇ ⋅ + ⋅ ∇ ⋅ ∇ ⋅ + ∇ ⋅ ⎡ ⋅ ∇ ⎤ + ∇ ⋅ ⎡ ⋅ ∇ ⎤⎣ ⎦ ⎣ ⎦ = − ⋅ ∇ ⋅ ∇ ⋅ + ∇ ⋅ ⎡ ⋅ ∇ + ⋅ ∇ ⎤⎣ ⎦ = + ww ww ww ww νν νν νν νν νν http://www.paper.edu.cn - 11 - ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) T T TT T1 1 1 T T TT T 2 2 2 u u u u u u λ λ λ λ λ λ ∇ ⋅ ⎡ ⋅ ∇ ⎤ + ∇ ⋅ ⎡ ⋅ ∇ ⎤ = ∇ ⋅ ∇⎣ ⎦ ⎣ ⎦ ∇ ⋅ ⎡ ⋅ ∇ ⎤ + ∇ ⋅ ⎡ ⋅ ∇ ⎤ = ∇ ⋅ ∇⎣ ⎦ ⎣ ⎦ νν νν ww ww 都是具有增强边缘作用的冲击滤波算子。同时，注意到当 0, 0σ ρ> > 取值相对较小时，下列 不等式 ( ) ( ) (