（计算数学专业论文）基于pde的变分水平集图像分割方法.pdf

独创性说明 作者郑重声明：本硕士学位论文是我个人在导师指导下进行的研究工 作及取得研究成果。尽我所知，除了文中特别加以标注和致谓t 的地方外， 论文中不包含其他人已经发表或撰写的研究成果，也不包含为获得大连理 工大学或者其他单位的学位或证书所使用过的材料。与我一同工作的同志 对本研究所做的贡献均已在论文中做了明确的说明并表示了谢意。 大连理工大学硕士研究生学位论文 大连理工大学学位论文版权使用授权书 本学位论文作者及指导教师完全了解“大连理工大学硕士、博士学位论文版权使用 规定”，同意大连理工大学保留并向国家有关部门或机构送交学位论文的复印件和电子 版，允许论文被查阅和借阅。本人授权大连理工大学可以将本学位论文的全部或部分内 容编入有关数据库进行检索，也可采用影印、缩印或扫描等复制手段保存和汇编学位论 文。 作者签名：堂塾堂 导师签名 大连理工大学硕士学位论文 1 绪论 数字图像处理，简称图像处理，即利用计算机对图像进行分析和处理，这一技术是随着 计算机技术的发展而逐渐发展起来的一个新的应用领域，汇聚了光学、电子学、数学、摄 影技术和计算机技术等学科的众多方面 1 】。图像处理中的许多问题，如图像分割，图像的 多尺度表示及图像重建等都产生了许多具有挑战性的数学问题。为便于计算机进行处理， 传统的图像处理方法都是先将图像离散化，得到一个二维图像矩阵，然后进行各种线性 运算。然而，随着图像处理在各个领域得到越来越广泛的应用，传统的图像处理方法的局 限性严重影响了图像处理的效果，难以达到人们所期望的目标。于是，人们开始尝试使用 更复杂的数学方法解决图像处理中更困难的问题。近年来，基于偏微分方程( p d e ) 的非线 性图像处理方法受到人们越来越多的关注f 2 ，3 。由于自然界的物体是连续的，而且图像本 身也是连续的，因此可以选择一定的函数空间，利用一个偏微分方程来描述物体成像的过 程。根据这一原理，便产生了基于偏微分方程的图像处理方法。 1 1 图像的定义 图像是利用各种观测系统以不同形式和手段观测客观世界而获得的，可以直接或间 接作用于人眼并进而产生视知觉的实体【4 。它是表达视觉信息的一种物理形式，是一种 高维的以人类视觉为最终归宿的信号。图像中的“图”就是物体反射或投射光的分布，“像” 是人的视觉系统接收图的信息而在人脑中形成的印象或认识。前者是客观存在的，而后者 是人的感觉，图像则是两者的结合。在一般意义下，一幅图像是另一个事物的一种表示，包 括了有关其所表示物体的描述信息。图像有各种各样的形式，其中一个重要的子集是由连 续函数和离散函数组成的抽象的数学图像，而后一种就是可以被计算机直接处理的数字 图像。 1 2 数字图像和数字图像处理 数字图像是指一个被采样和量化后的二维函数( 该二维函数由光学方法产生) ，采用 等间隔矩形网格进行采样，并对幅度进行等间隔量化，从而得到一个被量化的采样数值的 二维矩阵。图像的原本状态是连续的，处理之后的结果一般也要以连续的形式演绎，图像 的离散化表示只是为了便于利用计算机来实现算法 4 】。 所谓图像处理，是指对图像信息进行某种目的的加工修正以满足人的视觉心理或应 用需求的行为。从广义上分，图像处理可分为图像的模拟处理( 使用光学方法) 和图像的 数字处理( 使用电子学和数字方法) 。前者已有很长的发展历史，从简单的光学滤波到现 常先堂：基于p d e 的变分水平集图像分割方法 在的激光全息技术。后者通常称为数字图像处理，将图像信号转换成数字格式，并利用计 算机进行数学运算和分析，可以提高图像的实用性。数字图像处理技术的优点是处理精度 比较高，而且还可以通过改进处理软件来优化图像处理效果。自上世纪6 0 年代起，随着电 子技术和计算机技术的不断提高和普及，数字图像处理进入了一个飞速发展的时期。特别 是从上世纪9 0 年代开始，随着计算机多媒体技术和信息技术的蓬勃发展，突破了限制数字 图像处理发展的数据量庞大、处理速度相对较慢、硬件成本较高等的限制，数字图像处理 技术已逐渐深入到军事、商业、医学、艺术等领域甚至人们的日常生活中，并且在科学工 程领域也发挥了巨大的作用。目前，图像处理主要关心以下几个方面：图像预处理、图像 分割、图像边缘提取、形状建模、特征提取、目标识别、运动目标检测及图像可视化等。 1 3 图像分割及图像分割的意义 在对图像的研究和应用中，人们往往仅对图像中的某些部分感兴趣，这些部分通常称 为目标或前景( 其他部分称为背景) 。它们一般对应于图像中特定的、具有独特性质的区 域。为了辨识和分析目标，需要将它们提取分离出来，在此基础上才有可能对目标进一步 利用。图像分割就是指把图像分成各具特性的区域，并提取出感兴趣目标的技术和过程。 图像分割在不同的领域有时也用其他名称，如目标轮廓技术，阈值技术，目标检测技术，目 标识别技术等，这些技术本身或核心实际上也是图像分割技术【5 】。 图像分割是由图像处理到图像分析的关键步骤，在图像工程中占有重要位置。一方 面，它是目标表达的基础，对特征测量有重要的影响，另一方面，因为图像分割及其基于分 割的目标表达、特征提取和参数测量等将原始图像转化为更抽象、更紧凑的表达形式， 使得更高层的图像分析和理解成为可能。 1 4 图像分割的方法 1 4 1 图像分割的传统方法 传统的图像分割方法，一般的可分为以下三种 5 】： 基于阈值的分害0 方法 基于边缘检测的方法 基于区域提取的方法 其中，基于阂值的分割方法在本质上也是一种区域提取方法。由于图像种类的多样性和 复杂性，基于阈值的分割方法往往很难确定合适的阚值大小。而基于方向导数的边缘检 测方法对于噪声图像、边缘模糊图像或纹理图像的分割效果不理想，并且这种方法往往 需要一些预处理手段，如用高斯滤波去除图像中的噪声和预后处理，如进行边缘连接等。 基于区域提取的分割方法，如区域增长法，往往比较复杂。 大连理工大学硕士学位论文 1 4 2 基于p d e 的图像分割方法 为克服传统的图像处理方法的缺陷，近年来人们提出了应用p d e 进行图像处理，特 别是图像分割的方法。同传统的图像处理方法相比，p d e 方法具有以下突出的优点 6 】： p d e 方法是一种直接对图像进行分析的连续模型。在基于p d e 的图像处理 方法中，离散的滤波表现为连续的微分算子，因而使得网格的划分，局部非线性分析 易于实现。 由于p d e 方法可以直接处理一些图像中的梯度，曲率等几何特征，便于建立各种数 学模型灵活表述，而且可以充分利用数值分析和偏微分方程计算理论中的一些成 熟的结果，能获得较好的图像质量，并具有一定的稳定性。 应用p d e 方法，使得各种图像处理方法的合成十分自然。例如用于滤波合成的 实际应用闺像平滑与边缘保持，使得同时去噪和保持边缘的方案成为现实。 正是由于以上优点，p d e 方法已成为图像分析的一种重要方法。目前，基于p d e 的 图像处理已经发展为一个结合数学形态学、变分法、逼近论和仿射几何等数学方法 的完整的理论体系【7 。基于p d e 的图像分割是p d e 方法在图像处理领域中最具有 代表性的应用。 基于p d e 的图像分割方法总体上说可分为三类：是基于图像边缘信息的分割 方法，二是基于区域信息的分割方法，三是基于两者结合的方法。基于图像边缘梯度信 息的方法，对图像边缘具有较好的局域化作用，对于对比度较好的图像具有不错的分 割效果，但这种方法不适用于边缘模糊或噪声图像。基于区域的分割方法，主要利用图 像区域的统计信息和区域同质性的特征进行图像的全局分割。第三种方法同时考虑 了图像的边界和区域特征，是一种较理想的分割方法。一方面，它克服了仅仅依靠边界 梯度进行分割的缺陷，另一方面，在某些图像的分割应用中又能充分利用边界梯度信 息可以达到比仅利用区域信息进行分割更理想的效果。 在基于p d e 的图象分割方法中，基于变分法和水平集方法的活动轮廓模型集中 体现了p d e 图像分割方法的优越性，目前在p d e 图像分割中是一个研究的热点问题。 1 4 2 ，1 基于活动轮廓模型的图像分割方法 活动轮廓模型一般的可分为两种，一是基于变分方法的参数活动轮廓模型，二是 基于水平集方法的几何活动轮廓模型。相比于传统的图像分割方法，活动轮廓模型应 用于图像分割，具有以下显著的优点 8 ： 由于活动轮廓模型在连续状态下实现，所以最终得到的图像边缘表示可以达到 亚像素的精度，这一点对于医学图像分割和应用具有非常重要的意义。 通过约束被提取物体的边缘是光滑的，并且融入其它关于物体形状的先验信息， 活动轮廓模型对图像中的噪声和边缘间隙都具有较强的鲁棒性。 常先堂：基于p 睫的变分水平集图像分割方法 使用光滑的闭合曲线表示物体的边缘。避免了传统的图像分割方法中的预后处 理过程，如边缘连接等，而且对于物体的形状分析和识别都具有重要作用。 1 4 2 2 基于变分水平集方法的图像分割 基于变分方法的参数活动轮廓模型很难处理活动轮廓线拓扑结构的自适应变化， 而基于水平集方法的几何活动轮廓模型往往又不是能量极小化模型，故不满足能量极 小值原理。变分水平集方法综合了这两种方法的优点，使得分割模型既是一个能量极小 化模型，同时又可以处理演化曲线拓扑结构变化的自适应性。相比于传统的由纯粹的 p d e 驱动的几何活动轮廓模型，基于变分水平集方法的图像分割模型可以自然地融入 图像区域信息及先验形状信息，而且可以很好地处理一些传统的几何活动轮廓模型无 法解决的问题。 1 5 本文主要工作 经典的c v 模型【9 】是一种典型的基于图像区域信息的变分水平集分割模型，该模型 充分利用了图像的全局区域信息，在图像分割应用中取得了令人满意的结果。但是由于该 模型仅利用了区域信息，丽没有很好地使用图像边缘信息，使得在某些图像豹分割中，图 像边缘定位的精确性不够高，另一方面，由于该模型没有考虑到水平集函数自身具有的 内在性质，为保证数值解法的稳定性和模型演化的精确性，在某些应用中还需要周期性地 对水平集函数进行重新初始化。李纯明等【l o 】从变分的观点出发，针对水平集函数满足的 条件，提出了一种新的变分模型，完全去除了水平集函数的重新初始化步骤。但由于该模 型只考虑了图像的边界梯度信息，对边缘模糊图像，噪声图像和纹理图像不适用。 基于上述两个模型各自的优缺点，本文提出了一种新的基于边缘和区域信息的变分 水平集图像分割模型，一方面继承了上述两个模型的优点，另一方面又克服了各自的缺 点。同时推导出模型演化的梯度下降流方程，并证明了其隐格式迭代解法是无条件稳定的， 所以在保证数值解法稳定性的基础上，又可以适当地选取较大的迭代时间步长，以加速曲 线演化，减少迭代次数并提高了计算速度。 文章在第二章介绍了p d e 图像处理中两种常用的方法一交分法和水平集方法；第 三章对基于变分和水平集方法的活动轮廓模型在图像分割中的应用进行详细综述；第四 章是对本文工作的详细介绍，首先介绍了变分水平集方法，然后介绍了两种典型的变分 水平集图像分割模型，并分析了各自的优缺点，在此基础上提出了一种改进的模型，通过 不同的实际算例验证了本文方法的有效性。最后是对全文的总结和对未来工作的展望。 大连理工大学硕士学位论文 2 p d e 图像处理中两种常用的方法一变分法和水平集方法 2 。1 变分方法 变分方法是p d e 图像处理中一种常用的有效的数学方法。这一节将从泛函及求泛函 的极值出发，介绍变分法和梯度下降流方程，并以实例分析说明变分法在p d e 图像去噪的 t v 模型 1 1 中的应用。 2 1 1 变分法和泛函的一阶变分 2 1 1 1 泛函和变分法 ( 1 ) 泛函的定义：设兄r 分别是定义在同一数域p 上的赋范线性空间和实数空间，_ d 是x 的一个子集，若存在某种对应法则r 使得对任意x d ，有唯一确定的 y = r ( x ) = 疋y 与之对应，则r 称为泛函数，简称为泛函。这里，我们只考虑h i l b e r c 空 间中的泛函 1 3 】。 对泛函求极值的问题称为变分问题，使泛函取得极值的函数称为变分问题的解，也称 为极值函数，研究变分问题的学科称为变分法 1 2 】。 ( 2 ) h i l b e n 空间中的线性泛函；是指如下的一个映射z ：x 只，z 为h i l b e r t 空间，r 为 实数域，并满足如下条件：对任意的x ，y e z ，口，6 r 都有，( 甜+ 妙) = 口- ，( x ) + 6r f ( y ) 。若 存在常数c ，使得i f ( x ) i c 叫l ，则称线性泛函，( x ) 是有界的。对于一个给定的函数z x ， 若定义泛函，：0 ) = ，则易见它是一个有界线性泛函。对于h i l b e r c 空间中的有界线 性泛函有如下的定理。 定理1 ：设x 是一个h i l b e r t 空间，则对任意有界线性泛函z ：x r ，都存在一个函数z z ， 使得对所有的x 仨z ，z ( 对z ，x 成立。 2 1 1 2 泛函的导数( 一阶变分) 令f ：z 一丑是h i l b e n 空间z 中的一个泛函，设x x ，若下式 f ( 掣) ：陋坐竽型：丢即伽) b 、。 _ 0 五以 ” 的极限存在，则称它为泛函f 在x 处沿v 方向的方向导数。进一步，若f ( z ；v ) 是关于v 的有界线性泛函，则称，是g a t e a u x 可微的【1 3 】。根据定理1 ，由于f ( z ；v ) 是h i l b e r t 常先堂：基于p d e 的变分水平集图像分割方法 空司中的有界线性泛函，因此存在函数“肖，使得f 。( x ；v ) = 。称“是泛函， 的g a t e a u x 导数( 一阶变分) ，并记为“= f 。( x ) 。类似于函数求极值的必要条件，若， 是泛函f ( x ) 的极小解，则对所有的函数v z ，都有f + ( ，v ) = o 。特别地，f ( z ) = o 。 2 1 2p d e 图像处理中一类常用的泛函及其一阶变分 考虑泛函，( “) = p ( x ，“( x ) ，v “( x ) ) 出，其中工( x ，“，审“) 是关于变量x = ( 五，屯，矗) q ， “，v “的连续可微函数，且( x ，“，v “) r ( f 2 ) = 厂：j 1 厂( x ) 1 2 出 _ p ( “) 眈，其 中f ( “) 是f 的g a t ea _ 1 1 ) ( 导数( 一阶变分) 。由柯西许瓦兹不等式知，当v = 一f ( “) 时，上述积 分取得最小值，即泛函f 变化最快的方向是v = 一，( “) ，这个方向称为泛函f 的最速下降 方向t 称爱= 一f ( “) 为梯度流( 最速下降流) 方程，它描述了函数扰朝着泛函f 局部极 小值的变化。根据上述定理2 可知，对于能量泛函，( “) = 弘“( z ) ，v “( 砷) 出，其相应的梯度 流方程为： 詈= 喜毒喏( x ，扰，v 蝴一薷( 墨“，v “) ( 2 2 ) 下面，利用上述介绍的理论对p d e 图像处理中一个经典的用于图像去噪的总变差模型的 变分问题进行分析说明。 例1 ：总变差t v ( t o t a lv a r i a 士i o n ) 模型的变分描述为 1 1 】： 馏行f ( = 且甜( x ，力一，( 五) 2 矗q + j 1 v 甜( 五y ) l 矗q 其中，( x ，y ) 表示含噪声的图像，“( x ，y ) 为去噪处理后的图像，是权系数。 记w l = 鲁，w 2 = 等，令三( x ，”，v “) = 似一驴+ 段研+ 以，则 筹趔州) ，赞2 赢，奄2 赢 因此，根据( 2 1 ) 式，( “) 的g a t e 叽x 导数( 一阶变分) 为 f ( “) = 筹( 掣，v “) 一啬喏甜，v “) ) _ 2 ( 舻俨击( 静一号( 南 = 2 ( “一，) 一咖( 尚) 所以，求解泛函f ( “) 的极小值的梯度下降流方程为： 詈一( 咖枷嘲( 铲，) 常先堂：基于p 工) e 的变分水平集图像分割方法 2 2曲线演化理论和水平集方法 2 2 1 曲线演化的理论 由微分几何学的知识，我们知道描述瞄线几何特征的两个重要参数是曲线的单位法 矢丙和曲率七。其中，单位法矢丙描述曲线的方向，而曲率| 则描述曲线弯曲的程度。曲 线演化理论就是仅利用曲线的单位法矢和曲率等几何参数研究曲线随时间的变形，而这 些几何参数是与曲线的参数化方式无关的。 设平面闭合曲线c ( a f ) = ( 工( p ，f ) ，j ，( a r ) ) ，p 是任意参数化变量，是时间参数，并设 曲线的内向单位法矢为丙，曲率为七，则曲线沿其单位法矢方向的演化过程可用如下偏 微分方程表示簪= y ( c ) ( 2 。3 ) 其中，矿( c ) 是曲线演化的速度函数，决定曲线c 上每点的演化速度。如图2 1 中的箭头所示 为曲线的单位法矢丙，c 上各点的演化就是沿其单位法矢方向。e p s t e i n 和g a g e 已证明曲 线沿任意方向的运动总可以重新参数化为方程( 2 1 3 ) 的形式 1 4 1 。 图2 1 平面曲线演化 f 培2 1p l a r i 盯c u r v ee v 0 1 u t i o n 研究中最常用的曲线演化方式是曲率演化和常量演化。曲率演化可由如下的偏微分 方程描述 管= 口丘丙 ( 2 4 ) 其中，a 是正常数，j 是曲线的曲率。上述方程又称为“曲线变短流”( c u n ，es h o n e n i n g f l o w ) 。闭合曲线在它的作用下将会逐渐光滑地收缩变形，并使曲线的总长度减小，最终 使它收缩为一个圆点，如下图2 2 所示。 大连理工大学硕士学位论文 ( b )( c )( d ) 图2 2 曲线的曲率演化 f i g 2 2c u r ”e v 0 1 埘o nb 雎e do nc u r v d n j r ed e p e n d e ms p e e d 常量演化可由如下偏微分方程描述： 警= 哥 ( 2 - 5 ) 其中，k 是决定曲线演化速度的常数。图2 3 描述了曲线的常量演化过程，可以看到，常量 演化最终会导致曲线出现尖点，并可能出现拓扑结构的变化( 分裂或合并) 。同时，这个方 程又称为“面积减小( 增大) 流”，闭合曲线在这个方程的作用下将逐渐减小或扩大( 取 决于圪的符号) 它所包围的内部区域的面积，如图2 4 所示的效果a 唿甑唿8 。眨 ( a )( b ) ( c ) 图2 3 闭合曲线的常量演化( k = 1 ) f i 备2 ，3c o 璐t a n tf b ee v o l m i o no f c l 0 8 e d c u r v e s ( a ) ( b )( c ) 圈2 4 面积减小流( k = 1 ) f i g 2 4a r e ad e c r e a s i n gf l o w 下图2 5 为曲线基于曲率演化的“曲线变短流”和基于常量演化的“面积减小流”的对比。 9 常先堂：基于p d e 的变分水平集图像分割方法 ( a ) “曲线变短流”( b ) “面积减小流” 图2 5 “曲线变短流”与“面积减小流”的对比 f i g 2 5c o m p 撕s o n so f c u r v es h o r t e n i n g n o wa n da 碍a d e c r e a s i n g n o w 综上，我们得到了平面闭合曲线的两种常用演化方程，此时的演化曲线又称为活动 轮廓线。虽然我们可以使用参数化的方式表达曲线的动态演化，然而使用参数化的方式 计算曲线的曲率和法向矢量比较麻烦，并且需要处理演化曲线拓扑结构变化的情形，这 是一个相当复杂的处理过程，尤其是对三维情形。因此，需要有一种更加自然的处理曲线 的几何参数的计算方法，而水平集方法正好满足这个要求。 2 2 2 水平集方法 水平集方法( l e v e ls e tm e t h o d ) 【1 5 ，1 6 】最初由o s h e r 和s e t h i a n 提出，用于解决遵循热 力学方程下的火苗的外形变化过程。在水平集方法中，平面闭合曲线c 被隐含地表达为 三维连续函数曲面y ( x ，y ，f ) 的一个具有相同函数值的同值曲线，通常是 y = 0 ) ，称为零 水平集，妒( t 弘f ) 称为水平集函数。在处理曲线演化的问题时，按照一定的规律在坐标平 面上不断更新演化水平集函数，再求出零水平集所在的位置，即为曲线演化后的形状。 假设妒( x ，弘f ) ：r 2 o ，r ) 斗r 为l i p s c h i t z 连续的水平集函数，f 为时间参数，闭合曲 线族c ( b r ) ：o s p l 表示f 时刻与l c ，他弘f ) 对应的零水平集，即 c ( p ，o ) = ( 而力i 妒( 墨_ y ，o ) = o i c ( p ，r ) = ( 五力i y ( x ，y ，f ) = o 为使水平集函数y 在演化过程中，其零水平集所对应的平面闭合曲线矿( c ( f ) ，f ) = o ( z 6 ) 始终满足曲线演化的偏微分方程( 2 3 ) ，对方程( 2 6 ) 求全微分，得v y 等+ 警= o ( 2 7 ) 设s 是曲线c 的弧长参数，根据水平集函数的定义，y 沿着c 的切线方向的变化量为零， 即虬= o = 虬+ p ，儿= 。这样，v y 与闭合曲线c 的切线垂直，因此v i l f ，和c 的 大连理工大学硕士学位论文 法线同向。假设函数j ：f ，位于曲线c 内部的部分为负数，外部为正数，则零水平集c 的内向 单位法矢量为丙= 一斋，将它和( 2 3 ) 式代入( 2 7 ) 式，可得 詈一w 甲( c 污却”以c ) 尚卅( c ) l ( 2 8 ) 这个方程称为水平集方程 1 5 ，它是曲线演化方程( 2 3 ) 的欧氏表达。此时，曲线的曲率可 以用水平集函数表示为意= d 加崭) 。 下图2 6 可以帮助我们理解水平集方法的思想：假设一条圆曲线以依赖曲率的速度进行演 化，初始时刻的曲线如图2 6 ( a ) 所示，嵌入在曲面中的情形如图2 6 ( b ) y ( 五_ y ) = 0 ，到达另一 时刻曲线演化为图2 6 ( c ) ，此时嵌入在演化曲面中的情形如图2 6 ( d ) y ( 毛y ) = o ，从中可以 看到曲面y 也在不断地演化。 “，) 图2 6 演化的水平集及嵌入的曲面 f 嘻2 6e v 0 1 v i n gi e v e ls e tc u r v e 柚de 加血e d d e ds 删k e 常先堂：基于p d e 的变分水平集图像分割方法 一般的，水平集函数都是取为由初始曲线c n ( _ y ) 生成的符号距离函数( s d f ) 。设 矿( x ，弘o ) ，y ) r 2 是符号距离函数，则有y ( x ，y ，o ) = 如f ( x ，y ) 。其中d 衙( x ，y ) 表示 点( x ，y ) 到曲线c 0 的距离。对于任意闭合曲线，上述直接计算符号距离函数的计算量较 大。一般情况下，初始曲线都取为规则的闭合曲线，如圆或矩形等。如何快速准确地计算 任意闭合曲线的符号距离函数，对于提高水平集方法的效率和稳定性非常重要。s e t l l i a l l 等【1 7 】提出利用快速行进法计算符号距离，分两个步骤来完成，时间复杂度为o ( n l o g n ) 。 李俊等 1 8 】提出了一种基于源点映射的改进的扫描算法快速构造符号距离函数，整个计 算过程只需对图像网格点进行四次扫描，即可准确地计算出所有网格点的距离函数值， 计算复杂度仅为o ( n ) ，而且精度较高。 在水平集方法中，另一个值得注意的问题是速度函数的扩展问题。因为水平集方程 ( 2 8 ) 仅是从零水平集推导而来，一般的，速度函数只在零水平集上有定义，其他水平集上 未必有定义。因此，演化水平集函数时，需要将速度函数推广到所有的水平集，成为扩展 的速度场。许多扩展速度场的方法常常导致水平集函数不再保持为符号距离函数，从而 引起计算误差。如何选择合适的方法扩展速度场也是实现水平集方法时需要考虑的问 题。一种最简单的方法是直接计算每个网格点p 到初始曲线的最短距离的点，并以该点的 速度作为p 点的速度【1 9 】，但这种方法的计算量较大。a d a l s t e i n s s o n 等 2 0 】利用快速行进法 同时构造符号距离函数和扩展速度场，这种方法的最大优点是可以一次完成距离函数和 速度场的构造，但其计算复杂度仍然较高。 用水平集方法实现活动轮廓线的演化，实质上是把一个低维模型嵌入到一个高维空 间进行描述，虽然这种转化使得问题在形式上变得更加复杂，但在问题的求解上有许多 突出的优点： 只要速度函数y 是平滑的，则水平集函数妒( 五弘r ) 将始终保持为一个有效的函数， 而其零水平集即演化曲线c 的拓扑结构的变化能得到很自然的处理。 借助于偏微分方程中双曲守恒型的数值解法可以得到稳定的、高精度的有限差分格 式，并且可以得到唯一的、满足“熵守恒”条件的弱解。 水平集方法可以直接推广到求解高维曲面的演化运动。如利用三维闭合曲面的演化 可以对三维图像进行分割或对物体几何外形进行三维重建等。 2 2 3 水平集方法的数值计算 为实现曲线演化的水平集方法，必须采用适当的数值方法近似求解。由于水平集函 数在演化过程中始终保持为一个函数，因此可以用离散网格表达水平集函数矿( 五y ，f ) 。 2 一 设离散司格的间隔为| i z ，在，2 缸时刻网格结点( f ，) 处的水平集函数为瞄= y ( 访，舻，h 础) 这里& 为时间步长，则水平集函数妒的演化方程( 2 8 ) 可以离散化为： 警= w 啉i ( 2 9 ) 其中，w 表示拧出时刻扩展速度场位于网格点( f ，) 处的值。 在计算曲线演化方程中出现的法矢量和曲率时，需要用到水平集函数的一、二阶导 数，可以利用水平集函数的满足“熵守恒”条件的有限差分格式计算，以避免常量演化中 可能导致的水平集函数的奇异性问题。定义如下的六个有限差分算子： 叼= 警，叼= 警，班= 监铥 硝= 红蔷垃，叼= 纽，础= 警 则方程( 2 9 ) 的近似解可写成如下形式： 妒嚣1 = 嵋+ r d 【n a x ( k ，o ) v + + m i n ( k ，o ) v 一】 ( 2 1 0 ) 其中： v + = 【m 叫叼，o ) 2 + m i n ( 硝，o ) 2 + m 麟( 吲，o ) 2 + m i n ( 硝，o ) 2 】2 v 一= 【m i n ( d 才，o ) 2 + m “( p 孑，o ) 2 + m i n ( 日了，o ) 2 + m a x ( p 苫，o ) 2 】2 应用于图像分割，方程( 2 8 ) 中的速度函数通常可以写成如下形式：矿= ，+ k 。+ k 。， 其中= 是常最扩展速度，= 埘是曲率演化速度，= 矿开是水平对流速 度，哥= ( “，v ) 。这样，方程( 2 9 ) 的近似解可以写成下述形式： 矿嚣1 = 屹+ 址 【m a f ，o ) v + + m i n ( u ，0 ) v 一 m a x ( ，o ) 嘎j + m i n ( 1 易，o ) 叼 + m a x ( 吧，o ) 呸f + m i n ( 吗，o ) 叼】 一占 ( 功) 2 + ( j d 学) 2 弘 ( 2 1 1 ) 这种差分方法( 2 1 0 ) 或( 2 1 1 ) 称为逆向有限差分法( 迎风差分格式) 。 通过上述差分方程，可以利用迭代法不断地迭代更新水平集函数，然后利用轮廓检 测的方法，提取更新后的水平集函数的零水平集，即可得到演化后的活动轮廓线。在迭代 过程中，为了防止轮廓曲线发生位置漂移”，水平集函数应保持为符号距离函数，然而由 于采用数值解法，方程( 2 9 ) 的解并不总能保持为符号距离函数，因此需要周期性地对水 平集函数进行校正，即进行重新初始化的过程，以使它保持或接近符号距离函数，这是传 常先堂：基于p d e 的变分水平集图像分割方法 统水平集方法的一大缺陷，利用第四章中介绍的变分水平集方法，可以克服这一缺点，完 全去除了水平集函数的重新初始化过程，从而大大简化了计算过程。 2 24 水平集方法的快速计算方法 同传统的线性图像处理方法相比，水平集图像处理方法的一个最大的缺点是计算量 大。为了克服这一问题，人们提出了一些快速算法，常用的快速演化水平集的数值计算方 法有以下几种：第一种是由a d a l s t e i n s s o n 和s e 廿1 i a i l 提出的窄带水平集方法【2 1 】，第二种是 由s e t l l i a n 提出的快速行进法 1 7 】，第三种是由p a r a g i o s 等提出的h e 眦e s 算法 2 2 。 ( 1 ) 窄带水平集方法( n 删b 蛆dl e v e ls e tm e t l l o d ) 在传统的水平集方法中，由于每次都要对图像域上的每个网格点进行迭代计算，所以 计算量比较大。窄带水平集方法的主要思想就是在水平集曲面上，零水平集曲线的邻域内 定义一条具有一定宽度的窄带区域( n a r r o wb 髓d ) ( 如图2 7 ) ，每次只需更新计算窄带内 部区域中的网格点，通过不断更新窄带和零水平集进行求解。由于窄带中的点的数量相对 较少，所以这种方法可以大大减少计算量。 窄带水平集方法的具体算法如下： 初始化符号距离函数，依给定的窄带宽度和初始零水平集得到初始的窄带。 利用速度函数扩展的方法，计算窄带区域内网格点( f ，) 处的速度函数值。 利用迭代格式( 2 1 1 ) 在窄带区域内演化迭代得到妒。 使用轮廓线跟踪算法得到零水平集缈= 0 ，判断更新后的零水平集曲线与窄带 边缘的距离，若小于给定的闽值，则重新初始化符号距离函数和窄带，否则转。 由上述方法的计算过程可以看出。窄带水平集方法的计算效率与窄带宽度的选取有 很大关系。若窄带选取过窄。则重新初始化的次数就会增多，从而影响计算速度：若窄带 选取过宽，则每次迭代需要计算的点增多，也会影响计算效率。在实际计算中，往往根据具 体问题的不同和经验估计选取窄带宽度，这在一定程度上限制了窄带水平集方法的通用 性。 图2 7 窄带水平集方法的初始水平集曲面 f i g 2 7i 芏l i t i a ll e v e ls e ts u r f 如ei ni l 锄wb a l l dl e v e is e tm e m o d 大连理工大学硕士学位论文 ( 2 ) 快速行进法( f a s tm a r c h i n gm e t h o d ) 快速行进法是考虑曲线演化的某种特殊情形，即演化曲线的法向速度矿在演化过程 中始终保持符号不变，即总为正值或负值，此时可以把求解每个网格点的水平集函数的 问题等价地转化为求每个点的到达边界的时间r ( x ，y ) 。此时r ( 毛y ) 满足如下方程： l v r i 矿= 1( 2 1 2 ) 方程( 2 1 2 ) 是一个e i k o n a l ( 短时矩) 方程。利用迎风差分格式，它的稳定解可由如下方程得 到 一( ：篡：曲( 聊? 2 ， k ：上 ( 2 1 3 ) l + m a x ( d 了r ，o ) 2 + m i n ( d 孑r ，o ) 2 j 巧， 其中，叼r = 士生篙竽，端r = 丑参生分别是向前和向后有限差分。 快速行进法也是一种基于局部窄带区域的方法，与窄带水平集方法类似，先在图像域 内定义一条具有一定宽度的窄带，然后根据每个网格点距离初始轮廓线c 0 的距离，将网 格点分成三类： 激活点( a l i v e ”) ，曲线传播经过的点，到达时间r 已知。 活动点( “a c 廿v e ) ，在激活点周围的窄带区域中，瞳线传播即将经过的点。 远离点c f 嘣1 w a y ”) ，远离初始曲线的点，达到时间r 未知。 快速行进法的基本步骤如下： 1 1 初始化 将初始轮廓线c n 所在的网格点标记为“a l i v e ”，到达时间r = 0 。 考察每个a 1 i v e ”点的四邻点，如果有标记不是a l i v e ”的点，则标记为“a c t i v c ”， 并赋予到达时间为巧，= 彤，将所有标记为“a c t i v c ”的点放入一个以到达时间r 为序的二叉排序堆中。 剩余的既不是“a l i v e ”也不是“a c t i v e ”的点全部标为“f a r a w a y ”，到达时间为 巧j = 佃。 2 ) 行进过程 从排序堆中输出丁最小的标为“a c t i v e 的点j p ( f m 。矗。) ，将其标为a l i v e ”，并从排 序堆中删去r 值最小的点。 更新p ( f m 。厶；。) 的四邻点：将四邻点尸( k l ，矗，。) ，p ( f m m 十1 ，l t 。) ，p ( f m ”l 。一1 ) ， 常先堂：基于p d e 的变分水平集图像分割方法 p ( 。，l + i ) 中标记不是“a l i v e ，的点p ，根据( 2 1 3 ) 式计算其到达时间r 。若p 标为 “a c t i v e ”，则更新其到达时间，并调整在排序堆中的位置，若p 标为f a r a 、v a 咿，则将 其标记为“a c t i v e ”，并插入到排序堆中。 收敛检查：如果强。 o 是一个权系数，用于在向量扩散项与梯度符合程度之间进行折中。一般地，当图像 中噪声较多时，可适当地取个较大值。 相比于传统的s n a k e 模型，c f s n a k e 模型的优点有： 有效地扩大了外力的作用范围，对初始轮廓线的位置要求不太苛刻，如3 5 ( a ) 所示。 活动轮廓线可以收敛到物体的下凹区域( 如u 型谷的谷低) 如图3 5 ( b ) ( c ) 所示。 图3 5g v f s n a k e 的优点 f i g 3 5n ea d v a n t 8 9 e so f g v f s n a k e g v f 模型的缺点是g v f 外力的计算问题，通常采用显式迭代方法，但收敛速度较慢。 针对这一缺陷，林涛 3 6 利用多重网格的方法求解一个线性方程组，大大加快了g v f 外力的求解速度。另方面，g v f 的能量函数中是一个常数，这使得在图像边界附近 能量函数中扩散项的影响仍然较大，表现在最终曲线的演化结果上，就是曲线将会在图 像边界周围来回震荡。针对这一缺点，x u 等又提出了g g v f 模型 3 7 。g g v f 模型的 主要方法是将g v f 中的常数和i w f 分别替换为随图像边界梯度l v 厂i 的大小而变化的 函数，即令g ( 1 w 1 ) = e x p ( 一1 w i 茧) ， ( i 可1 ) = l g ( 1 吖 ) ，k 为自定义的常量，这两个函数用 来在扩散场的光滑度和梯度符合程度之间进行折中。求解g g v f 的过程如下 粤 = g ( 1 1 盯i ) v 2 矿一 ( i 可1 ) ( 矿一1 9 力 ( 3 7 ) 初值条件为矿( x ，y ，o ) = 盯。o g v f 模型的优点是，对于下凹程度较大的物体边缘也具有 较好的收敛效果，如下图3 6 ( a ) ( b ) 所示。但g g v f 模型也有一个缺点，即它并不是从变 分原理推导出的，故不满足能量极小值原理，而且显式迭代数值解法的收敛速度也很慢。 常先堂：基于p d e 的变分水平集图像分割方法 图3 6 g g v f s n a l c e 模型 f 培3 6g g v f - s n a l 【em o d e l 林涛 3 6 在总结了g 、，f 和g g v f 模型各自的优缺点后，提出为进一步提高计算速 度并且保证隐式迭代格式的收敛性，采用普通g v f 的迭代方法，将g v f 方法的迭代过 程转化为隐格式求解。并且，为了克服g v f 在图像边缘附近不准确的缺点，将得到的 g v f 能量模型的收敛解矿( x ，y ) 与图像边界梯度场逐点进行加权求和，即令 矿= g ( i v ，i ) 矿+ ( i v ，i ) 可这样，既保证了在图像边界附近与边界梯度的一致性，又可将图 像原先作用范围狭小的边界梯度的影响力扩大到整个图像范围。 ( 5 ) 动态距离力 基本思想是通过计算变形曲线上每点处的符号距离求出一个外力。每点处的符号 距离等于沿该点法线方向到距其最近的边缘点的距离。这个符号距离值每当轮廓线演 化一次就需重新计算。此时，要设计一项标准判断哪些像素点是图像边缘上的点。同 时，还需要设定一个阈值作为最大搜索距离，这样设计的外力称为动态距离力，它也 能将活动轮廓线从距离图像边缘较远处拉向图像的边界。这种外力的缺点也很显然， 进行一维搜索( 沿法线方向) 较耗时，且每当曲线演化移动一次后，都要重新计算。 3 1 3 3 参数s n a k e 模型中轮廓线的表达方式 对于参数化的s n a k e 模型，从模型曲线的表达方式上分，一般的可分为两种，即基 于离散点的表达方式和参数化的表达方式传统的s n a l 【e 模型就是一种典型的基于离散 点的表达方式，即活动轮廓曲线由离散的有序点列组成，活动轮廓线的演化就是对这些 离散点列进行演化。但这种表达方式有许多缺点： 由于需要同时优化大量系数，导致曲线的收敛速度非常缓慢 使用不连续的离散点列描述曲线，增加了计算复杂度 不易确定各光滑约束项的对应系数 在噪声影响较大时，离散曲线的高阶导数计算可能不准确 大连理工大学硕士学位论文 为解决上述问题，人们提出了活动轮廓线的各种参数化表达形式，目的都是为了使s n a k e 模型更加稳定和更快地收敛。其中最有代表性的是b s n a k e 模型 2 6 ，2 7 】。 b s n a k e 模型是一种将曲线的参数化b 样条表示和s n a k e 模型相结合的活动轮廓模 型。用参数化的b 样条曲线表达轮廓线有许多优点，如局部可控制性，光滑性等。b s n a k e 方法允许对变形曲线的局部控制，具有紧凑的表达，而且所需优化的参数个数较少，光 滑性由模型自身的属性直接保证，从而不需要特别的光滑项约束，可通过对个别控制点 的局部调整，进行s n a k e 模型的局部调节。 在实际应用中，具有极小曲率性质的三次b 一样条曲线是最常用的曲线表达方式。令 s 。( 功为节点间距为 的三次均匀b 一样条曲线，即s 。( x ) = g ( 女) 卢3 ( 三一女) ，其中矿( 曲表 示三次均匀曰样条基函数，g ( 七) 表示样条曲线的控制点组。由于在应用b - s n a k e 方法 时，演化曲线的光滑性由b 一样条曲线自身的属性保证，因此传统的参数化s n a k e 模型中 的内部能量项在b s n a k e 模型中可以省略，这样b s n a k e 模型中只有外部能量项。令 g ，) = 、f 馔+ ，) 2 + 鬻厂) 2 ，其中厂为图像函数，表示图像的光滑核函数。这样， b s n a k e 模型的外部能量项可写为e ( c ( 女) ) = 一g ( j ( m ，其中j ( f ) 表示b 一样条曲线上的 第f 个离散采样点，c ( d 表示样条曲线的第个控制点。为求得样条曲线的控制点组 f c ( ) 的演化方程，可使用如下的梯度下降法： 帮= 篙1 鬻，鬻= 等，嚣= 等c 黼 a c ( 女) 盘a c ( t ) 4a c ( 女)( f )a c 位)船( f ) 。” 7 下图3 7 所示为基于b s n a k e 模型的医学c t 图像分割的一个例子。 ( a ) 初始b 一样条曲线 ( b ) 分割曲线 图3 7 基于b s n a k e 的医学图象分割 f i g3 7m e d j c a ii m a g es e g m e t l t a t i o nb 如e d o nb s n a ；c e 2 7 常先堂：基于p d e 的变分水平集图像分害4 方法 3 2 几何活动轮廓分割模型 几何活动轮廓模型较好地克服了参数化s n a k e 模型的许多缺点，如可以处理演化曲 线拓扑结构的变化，对初始位置不敏感。具有稳定、唯一的数值解等。几何活动轮廓模 型已引起人们越来越多的关注，并已在图象处理和计算机视觉领域得到广泛应用。 3 2 1 几何活动轮廓分割模型的简单分类 ( 1 ) 基于曲线演化方法的几何活动轮廓分割模型 在图像分割的过程中，当我们给定了一条初始轮廓曲线，则可以利用图像数据的某 种属性来驱动轮廓线向图像中的目标边界演化，这样即得到了基于曲线演化方法的几何 活动轮廓图像分割模型。由第二章介绍的曲线演化理论知识知，曲线c ( p f ) 沿法线方向的 变形可用如下的p d e 描述：丝弘= 位七+ ) 帚。其中七为曲线的曲率，v o ，口是常数。这 里，口+ 相当于曲线沿法线方向的演化速度。由于演化曲线要在图像边界处停止，即演化 速度在边界处取极小值，因此可以在这个速度项的前面乘上一个与图像梯度相关的属性 函数g ( 1 w i ) = 1 干阿0 砑，这样便得到由c a s e l l e s 3 8 】和m a l l a d i 【1 9 】各自独立提出的基于曲 线演化方法的图像分割模型 筹= g d v 加 女+ v o ) ( 3 8 ) 这时，曲线的演化速度g a v 巾 _ j + v o ) 在图像边界处取极小值。该模型一般的称为几何 活动轮廓模型。利用水平集方法，将其表示成水平集函数的演化方程为 等= g ( 删) ( 口t + v o ) l v 训= a g ( 删) 咖搦，i v 卅v o g ( i v 巾i v 纠 ( 3 9 ) 下图3