第四章偏微分程数值求解以及运算结果展示.doc

第四章 偏微分程数值求解以及运算结果展示4.1偏微分方程数值求解和水平集方法用于处理图像问题的偏微分方程比较复杂，往往无法求出一个解析解,通常是求数值逼近解。因此，理论上说，偏微分方程数值求解的一般方法都可以拿到图像偏微分方程的求解中来使用。从原则上说，由于图像偏微分方程是利用微分几何的方法来描述一幅图像，以及一个图像处理问题。而图像信息主要蕴含在像素灰度级的跳变（图像边缘信息），以及表变得几何分布上。因此，图像偏微分方程往往具有其复杂性。而图像偏微分方程的求解，往往需要复杂的技巧性。常用的偏微分方程数值求解方法有:有限差分法,有限单元法，边界元法等等，其中，对图像偏微分方程求解主要使用有限差分法,迎风格式的有限差分法在图像偏微分方程求解中广泛地使用。目前在偏微分图像处理中，通常辅以使用Osher和Sethian提出的水平集方法（Level Set Method）进行数值求解。4。1.1图像偏微分方程数值求解水平集方法（Level Set Method）最初由Osher和Sethian提出，其主要思想是将曲线、曲面和图像演化表示为更高维数的超平面水平集，其演化速度是该曲线或曲面的局部曲率的函数。因此，水平集方法可以看作是使用欧拉方法求解隐性偏微分方程的一种具体实现方式，具体地说,以计算的复杂度换取对拓扑变化的适应性,能处理比较困难的曲线或曲面演化过程拓扑结构变化问题。水平集法一般需要有初始的轮廓线或轮廓面和速度图像。其中初始轮廓线或轮廓面一般由用户手工勾划出或输入种子点代替。速度图像可以根据图像的特征如曲率、平均曲率或梯度而定。水平集演化产生的图像是一个时间的函数图像,称为时间距离图（TimeDistance Plane）。即对每个象素，轮廓前沿到达此象素所经历的最短时间。通过指定一个停止阈值，当前沿进化迭代到一定数目时停止,这就产生了以演化时间为参数的等值面序列,直至完成图像处理任务。 初始化水平集方法时，需要区分两维闭合曲线的内外部，以构造符号距离函数(Signed Distance Function),而如何区分任意形状闭合曲线（闭合曲面）的内外部是一个比较烦琐的问题.另外，水平集方法迭代过程中,为了计算稳定，往往要间隔地对水平集函数重新初始化（Re-initialization），而常规的直接计算点与曲线距离的初始化方法5计算时间较长，大大增加了计算量.因此，如何有效而快速地进行水平集函数的重新初始化，对提高水平集方法的计算效率很有意义.水平集方法进行数值求解,在偏微分方程图像处理领域已经得到了广泛的应用。比如， Alvarez9 以及Malladi，Sethian5 等用水平集去除图像噪声；Parogios以水平集方法进行纹理分割以及运动目标分割7,8； Bertalmio 等将水平集方法应用于图像变形和破损图像修复(Inpainting)中10，11;Mansouri 将水平集运用于运动目标跟踪邻域12；Fuaeras和Keriven13以水平集方法解决立体匹配问题。在计算过程中，所有网格点上的函数值都有更新的计算.因此计算量往往是很庞大的，尤其是对于计算用于图像处理的偏微分方程时，更是如此.复杂的计算方法所换来的优势是，曲线的拓扑结构变化处于整个算法的控制之下.因为不论曲线如何变化，它始终处于平面或空间网格中，而演化后的轮廓跟踪步骤可以得到平面上任何位置的曲线。因此，水平集方法求解也大大促进了曲率方程曲线演化模型的迅速发展66。4。1.2水平集方法图像微分方程作的求解,通常使用的必要工具是图像的水平线(Level Set Contours），即，图像的水平集各自对应的边界线。在这个框架下，一幅图像就由一系列的零水平集函数表示。图像处理的过程就表现为相应的图像的水平线的移动/演化过程。因此，在流体力学中表达各种曲率流的方程在图像处理中得到了广泛的应用。 水平集方法中,图像偏微分方程可以用一个普遍的表达式来表示。一幅图像可以由一系列的零水平集函数表示。图像处理的过程表现为水平集函数组所约束的曲线（称为水平线）的演化过程。我们可以把以上所述的图像偏微分方程用一个普遍的表达式来表示。设：表示一幅灰度图像,而是灰度值。引入迭代次数参量（“伪时间”参量),图像的变化可以用偏微分演化方程：是随时间变化的图像，而是初始条件，是表示所使用的算子.偏微分方程的解给出了伪时刻时（即迭代次数为）图像的状态。对于矢量值图像，可以建立相应的偏微分方程组。 在这个框架下表示曲线(轮廓线,水平线）的演化方程。对于曲线是的函数.曲线的演化视为在法线方向上的变形，而变形速度与相应曲线或曲面的曲率有关，方程如下 =N 其中：是曲率， N是曲线的法矢量，是曲线演化算子。有时，方程中可能也同时出现切向速度，但这一速度分量不影响曲线几何形状的变化。这一方法可以方便地推广到曲面演化模型的情况下。 水平集方法（Level Set Method) 最初由Osher和Sethian提出28，其主要思想是将曲线、曲面和图像演化表示为更高维数的超平面水平集，其演化速度是该曲线或曲面的局部曲率的函数。因此，水平集方法可以看作是使用欧拉方法求解隐性偏微分方程的一种具体实现方式，具体地说,以计算的复杂度换取对拓扑变化的适应性,能处理比较困难的曲线或曲面演化过程拓扑结构变化问题。在计算过程中，所有网格点上的函数值都有更新的计算。复杂的计算方法所换来的优势是,曲线的拓扑结构变化处于整个算法的控制之下.因为不论曲线如何变化，它始终处于平面或空间网格中，而演化后的轮廓跟踪步骤可以得到平面上任何位置的曲线。Osher-Sethain提出的水平集方法有效地解决了在活动围线演化过程中曲线拓扑结构变化的问题。因此，水平集的方法在包括图像处理在内的诸多领域内，迅速地得到了广泛应用.但是传统的水平集方法对图像只能作两个分片区域的分割，而且在数值计算中常常需要引入对零水平集函数的重新初始化过程,造成计算的繁琐。针对这些弊端，适合图像分割的新的水平集方法不断被提出。典型的技术是VeseChan多相水平集方法和P-D测地活动区域组,其他的水平集方法还有1，4-6，23,40。4。2运算结果展示以及对于水平集方法求解技巧的讨论4。2.1运算结果展示本论文工作主要跟踪了偏微分方程图像处理方面的一些感兴趣的问题，尤其对偏微分方程图像分割和图像修补问题的研究方法和研究动态有了较为清晰的认识。本工作主要完成了两个方面的图像修补工作。本工作的第一个方面是利用ROF图像复原模型和BSBC图像修补模型，在图像子空间中进行修补的算法5，48。工作包括如下三个步骤： 改进R-O-F图像复原模型，从而使得原模型中被摒除的部分成为输入图像的纹理子图像。 基于改进了的RO-F图像复原模型进行图像分解,得到输入图像的结构子图像和纹理子图像。 在图像分解子空间中进行图像修补。其中对结构子图像做和纹理子图像。本实验中所用的测试图像有两个。一个是使用MATLAB合成的简单的几何测试图像，图像中也包含了简单的纹理信息。另一幅是自然图像白塔。所使用的彩色图像首先被分解为R，G,B，子图像，然后分别通过本修补程序处理；最后重新合成为彩色图像。几何测试图像的测试图像结果如下。 图41。1 原始图像, 红色圆形有缺损fig.41。1 input image 图41.2 中间结果，图像缺损部分正在被填补 fig。4-1.2 intermediate result 图41。3 输出图像，红色圆形缺损部分基本被填补fig。4-1.3 inpainted image自然图像白塔的测试图像结果如下。 图4-2。1 原始图像,白塔将被剔除 fig。4-2。1 original image, with tower to be removed 图42.2a 剪切原始图像，以节省计算量fig。4-2.2a tailored image to reduce calculation 图42。2b标记图像填补区域 fig。4-2.2b inpainting mask 图4-2。3 a,b，c R-，G，B子空间图像,分别作为本修补程序的输入图像fig。42。3 a，b，c R-，G-，Bsubspace of tailored image as input image of the algorithm 图42.4 a，b,c R,G，B子空间图像的中间结果fig。42.4 a,b，c intermediate result of R-,G,B-subspace of tailored image 图4-2。5a 输出图像图42.5b输出图像7fig。42。5a result image fig。4-2。5b result image of 7本工作的第二个方面是再现了EBI最佳相似复制法24。这是一个非基于偏微分方程的图像修补方法,其主要算法包括如下三个步骤： 对待填补区域边界上的像素点P进行优先级排序。其中，邻近的已知像素信息越丰富的点更具有优先权。 在图像的已知区域内寻找与待填补区域边界上的像素点P相对应的最佳相似点Q，并将相应的像素信息向待填补区域内部推进。 依次反复推进，直至修补完毕.过程中，有必要进行人工干预以保证最佳相似点的选取和像素点优先级排序没有出现明显的差错。实验中的测试图像结果如下。图4-3.1 原始图像,图像前景人物将被剔除fig.4-3。1 input image，woman to be removed 图4-3.2 输出图像，人物用图像背景成功地填补fig.4-3.2 image inpainted在附录中，列出了我用MATLAB语言实现的两个图像修补算法的源程序片段:BSBC在图像子空间中修补5，和一个非基于偏微分方程的图像修补方法，EBI最佳相似复制法24。4.2.2对于水平集方法求解技巧的讨论我的工作主要是试图理解偏微分方程求解的基本方法，并通过跟踪现有的方法和实现其提出的水平集方法辅助求解的基本策略，试图理解水平集方法在偏微分方程图像处理中的使用技巧.附录中的BSBC-偏微分方程图像处理程序，是利用平均曲率移动（Mean Curvature Motion)，或称为平均曲率流(Mean Curvature Flow）的模型，使用的是迎风格式（Upwind Scheme）的有限差分法。平均曲率移动是一个数学性质很好的偏微分方程,因而，比较容易用简单的有限差分法实现。 重新初始化（Reinitialization ）和粘性延伸（Velocity Extension），是现有水平集（Level Set)方法中两个重要的技巧.然而，这些都是为了保证曲线演化（Evolution）的稳定和可用的结果，人们不得不采取的补救措施（Remedy）。当然，绝大多数水平集的文章都不说这是缺陷，相反，这些补救措施的实现方法都成为一个重要的课题。比如说，单单关于如何做重新初始化这方面的文章就有很多。 所谓重新初始化（Reinitialization，顾名思义就是反复初始化。在迭代之前就要把水平集函数（Level Set Function）初始化为一个数学性质很好的函数.通常是利用符号距离函数（SDF， Signed Distance Function),然后每迭代几步就要停下来，修改这个函数，把它改成一个新的符号距离函数；再接着迭代，周而复始.这个每迭代几步就停下来修改水平集函数的过程就是重新初始化。 在修改水平集函数的同时，还往往要修改或重建速度函数（Speed Function）(就是所求解的水平集函数表示下的偏微分方程的系数)。这个过程叫粘性延伸（Velocity Extension）。 每做一次重新初始化，就需要重新计算符号距离函数（SDF， Signed Distance Function）,以确定各个像素点演化的优先顺序。符号距离函数，其符号是指图像中的一个像素点与零水平曲线的相对位置，由Heaviside函数确定；其大小是指该像素点到一个闭合曲线零水平曲线上的最短距离，这个问题比较麻烦(尽管有些文章也提到了一些算法). 如果是指给定封闭曲线，算相应的符号距离函数，还是很直接的.可是在对水平集函数进行重新初始化的时候,即每迭代几步，停止迭代，计算所得到的水平集函数的零水平集曲线（Zero Level Contours）的符号距离函数,这就不容易了.现在主要有三种方法: 1。 确定水平集函数（以下记为u)的零水平集曲线（这可以用Matlab函数contour得到），然后算这些零水平集曲线的符号距离函数。这种方法的思路很直接，但是在实际应用中，如果u的零水平集函数（Zero Level Set）由几个零水平集曲线构成，而且出现大零水平集曲线里面还有小零水平集曲线的情况，这时候在判断哪些点是这些零水平集曲线的内部外部(即确定SDF的符号)会相当麻烦。 2。 用偏微分方程（PDE，Partial Differential Equation）的方法，即解一个基于时间（time dependent）的偏微分方程，其平稳稳定解（state steady solution）就是一个符号距离函数。 3。 用Sethian最先提出的快速步进法,或称为快速遍历法（Fast Marching），得到符号距离函数50。本文工作使用了方法3的思路源距离扫描算法,四次扫描就可以完成一次距离函数的重新初始化，由于在整个图像定义域内更新水平集函数,所以有限几次迭代就可以得到收敛解66. 另一方面，重新初始化和粘性延伸是现有水平集方法中两个重要的数值求解补救措施。它们不仅增加了计算量,而且它们本身也没能克服全部的问题，对最后结果有着很复杂的影响,甚至是负面的影响.于是有的研究人员索性试图摆脱这种被动地采用补救措施的思路，通过修正偏微分方程模型等途径，实现比较简单的求解.试图在不影响演化的稳定性和最后结果的精确度的前提下，可以对时间步长的限制可以选得很大，从而加快演化。本文工作没有涉及这方面的内容。4。3本章总结本章首先阐述了对偏微分方程进行数值求解的一般方法,然后对于目前流行的利用水平集方法对应用于图像处理的偏微分方程进行数值求解的方法进行深入的探讨。 本章的最后是在本文工作中所积累起来的一些数值运算技巧和方法的总结和讨论，以及我实现出来的几个算法在基本测试图像上的运行结果。对应的MATLAB源程序在附录部分。本章没有述及一些研究人员所做的不需要引入数字图像的零水平集函数而进行直接求解的新的尝试.