（应用数学专业论文）level+set方程的求解与应用.pdf

浙江人学颇| 学位论文l e v e ls e t 方程的求斛。引训可 摘要 摘要 l e v e ls e t 方法将运动界面表示为高维场函数的零等值面，通过l e v e ls e t 方程 驱动场函数演化，从而跟踪运动界面的演化，自然而鲁棒地解决了界面演化中拓 扑结构改变的问题，在计算机图形学、计算物理、图像处理、计算机视觉、化学、 控制理论等众多领域得到了广泛应用。 本文以l e v e ls e t 方程的求解与应用为中心，简要回顾t l e v e ls e t 模型的理论 背景和经典的求解算法及其应用现状，同时对当前的研究热点进行了归纳总结， 讨论t l e v e ls e t 模型的优点和存在的问题，进而介绍我们在l e v e ls e t 方程求解以 及应用方面的工作。 求解动态l e v e ls e t 方程非常耗时，均匀、各向同性的方法简单易用，但不够 灵活，往往不能满足实际应用的需求。本文在前人工作的基础上，针对三维空间 中的三角网格曲面流形演化的特点，提出了基于窄带的自适应l e v e ls e t 方法。自 适应方法构建粗网格满足界面演化的整体需求，同时估算粗网格点的曲率值，使 用快速扩散法聚类高曲率点，构建细网格捕捉界面演化的细节。晟后我们将自适 应方法应用予几何模型的m o r p h i n g 动面。实验结果与误差分析表明，自适应方法 能有效减少计算量，达到更好的界面跟踪效果。 结合网格方法与矢量方法各自的优点，本文提出了混合等距线算法。算法求 解静态l e v e ls e t 方程生成距离场，从中抽取等距线的拓扑结构信息，通过顶点调 整、奇异处理、特征识别等步骤最终将等距线表示为基本图元的组合。基于流形 上的静态l e v e ls e t 方程求解，本文还实现了带边界三角网格曲面上的广义等距线 和测地线算法，进而生成参数网格计算原始曲面的平面参数坐标，达到了较好的 参数化结果。 关键词：三角网格，l e v e ls e t ，窄带，白适应，聚类，演化，等距线，图元，测 地线，参数化 浙江人学 ! i 。学位论文k v e ls e t 方程的求解j 成用a b s t r a c t a b s t r a c t t h el e v e ls e tm o d e lr e p r e s e n t st h ei n t e r f a c ea st h e0 一l e v e ls e to ft h ee m b e d d e d f i e l df u n c t i o na n dt h ee v o l u t i o ni sd r i v e nb yt h el e v e ls e te q u a t i o n t h i se l e g a n t p e r s p e c t i v ep r o c e s s e st h et o p o l o g i c a lc h a n g eo fi n t e r f a c ee v o l u t i o nn a t u r a l l ya n d r o b u s t l y i th a sb e e nw i d e l ya p p l i e dt om a n yf i e l d si n c l u d i n gc o m p u t a t i o n a lp h y s i c s ， c o m p u t e rg r a p h i c s ，i m a g ep r o c e s s i n g ，c o m p u t e rv i s i o n ，c h e m i s t r y , c o n t r o lt h e o r y , e l c t h et h e s i sa d d r e s s e st h et o p i c sa b o u tt h ee f f i c i e n ts o l u t i o no ft h el e v e ls e t e q u a t i o na n di t sa p p l i c a t i o n s t h em a t h e m a t i cb a c k g r o u n da n ds o m et r a d i t i o n a l s o l u t i o nm e t h o d sa r eb r i e f l yi n t r o d u c e d t h e nt h es t a t eo fa r t so nl e v e ls e tm e t h o d sa r e s u n a r n a r i z e da n da n a l y z e d ，m e a n w h i l eo u rt e c h n i c a lc o n t r i b u t i o n so nl e v e ls e tm e t h o d a r ed e s c r i b e d n u m e r i c a l l ys o l v i n gl e v e ls e te q u a t i o ni sa t i m ec o n s u m i n gp r o c e d u r e u n i f o r m a n di s o t r o p i c a p p r o a c h i s s i m p l eb u ts o m e w h a ti n e f f i c i e n t m a n yr e s e a r c h e s a r e c o n d u c t e do na c c e l e r a t i n gt h en u m e r i c a ls o l u t i o n t oe f f i c i e n t l yu a c et h e2 - d t r i a n g u l a rm e s hm a n i f o l di n3 - ds p a c e ，an o v e la d a p t i v el e v e ls e ts o l u t i o ni sp r o p o s e d w h i c hi sb a s e do nn a l t o wb a n dm e t h o d i tb u i l d sc o a r s er e s o l u t i o ng r i df o rg l o b a l m e s he v o l u t i o n ，t h e nh i g hc u r v a t u r eg r i dp o i n t sa r ed e t e c t e da n dc l u s t e r e du s i n gf a s t d i f f u s i o nm e t h o d t h ec o r r e s p o n d i n gf i n er e s o l u t i o ng t i d sa r eb u i l tt oc a p t u r et h e d e t a i l so nt h ei n t e r f a c e t h ea d a p t i v em e t h o di sa p p l i e dt ot h em e s hm o r p h i n g a n i m a t i o n i m p l e m e n t a t i o nr e s u l t sa n de r r o ra n a l y s i ss h o wt h a tt h ep r o p o s e da d a p t i v e m e t h o dh a sa d v a n t a g e so nc o m p u t a t i o n a le f f i c i e n c ya n da c c u r a c y a h y b r i dc u r v eo f f s e tm e t h o di sp r o p o s e db yc o m b i n i n gt h eg r i da p p r o a c ha n d v e c t o ra p p r o a c h t h et o p o l o g i c a li n f o r m a t i o no ft h eo f f s e tc u r v e sc a nb ee x t r a c t e d f i o mt h es o l u t i o no ft h es t a t i cl e v e ls e le q u a t i o n t h ei s o l e v e ls e tc u r v e sa l ef i n a l l y r e p r e s e n t e da sa s e to ft h eb a s i cg r a p he l e m e n t st h r o u g hv e r t e xa d j u s t m e n t ，s i n g u l a r i t y p r o c e s s i n ga n df e a t u r ed e t e c t i o np r o c e d u r e b yc a l c u l a t i n gt h e o f f s e tc u r v e sa n d g e o d e s i cc u r v e so nt h em e s h ，ad e s i r a b l ep a r a m e t e r i z a t i o no ft h em e s hw i t hb o u n d a l y 浙江人学倾i ：学位论文l e v e ls e t 方程的求解j 埘用a b s t r a c t i so b t a i n e d k e yw o r d s ：t r i a n g u l a rm e s h ，l e v e ls e t ，b a r r o wb a n d ，a d a p t i v e ，c l u s t e r , e v o l u t i o n ， o f f s e tc u r v e ，g r a p he l e m e n t ，g e o d e s i c ，p a r a m e t e r i z a t i o n i i l 浙江大学颁i ：学位论文k v e ls e t 方程的求解j 应用 第一节绪论 第一章绪论 运动界面的传播演化出现在我们生活的每一个角落，潮汐涨落，云雨升腾， 火焰燃烧等自然现象都是生动的实例。如何自然高效地跟踪运动界面演化，真实 再现各种自然场景是图形学领域一个重要的研究课题。1 9 8 8 年o s h e r 与s e t h i a n 提出了l e v e ls e t 模型，通过偏微分方程描述界面演化的动态过程，自然而鲁棒 地解决了界面演化中拓扑结构改变的难题，为后续研究工作奠定了重要的理论基 础。九十年代以来，许多学者相继加入了l e v e ls e t 方法的研究队伍，使得l e v e l s e t 方法被广泛应用于计算机图形学、计算物理、图像处理、计算机视觉、化学、 控制理论等众多领域。s i g g r a p h 会议分别于2 0 0 2 年和2 0 0 4 年出了两个专题研讨 l e v e ls e t 方法。如今，l e v e ls e t 方法正处于飞速发展的阶段，每年都涌现出大 量新颖的工作。 1 1 运动界面跟踪方法 研究运动界面的演化有许多种方法，通常可以归结为l a g r a n g i a n 观点和 e u l e f i a n 观点两种。l a g r a n g i a n 观点关注运动界面上每个质点的运动，而e u l e r i a n 观点则研究既定空间区域中界面的运动状况。 基于l a g r a n g i a n 观点的浮标法( m a r k e rp a r t i c l e ) 将运动界面离散成为采样点的 集合，跟踪采样点在运动规律作用下到达的空间位置，通过重建算法将采样点连 接成线段( 2 d ) 或三角片曲面( 3 d ) 。浮标法简单直观，但存在严重缺陷，界面演化 多个时间步长后，采样点的空间位景分布变得非常不均匀，无法提供足够的信息 重建出运动界面的空间位置；同时很难描述奇异点处的运动规律；更重要的是， 当界面发生分裂或合并等拓扑改变的情形，则很难正确地重构界面，算法也不容 易推广到高维空间。 基于e u l e r i a n 观点的v o l u m eo ff l u i d 方法( 简称为v o f 方法) 将研究区域划 分为离散的网格单元，每个单元按照包含材料的多少被赋予相应的密度值，位于 界面内部区域的单元密度值设为1 ，外部区域的单元密度值设为0 ，界面经过的单 浙江人学倾i 。学位论史k v e ls e t 方程的求斛j 晦用 第一章绪论 元密度值设为界面包围的体积与单元体积的比值。演化过程中根据界面的运动方 程修正各个单元的密度值，晟后抽取等值面得到运动界面的空间位置。v o f 方法 可以很好地处理拓扑结构改变、界面自相交等问题，但基于网格单元的v o f 方法 还存在较大的局限性，比如估算曲面的法向、曲率等内在几何量的误差较大；往 往需要将区域划分成大量的网格单元才能获得较高的计算精度；对奇异点的处理 方法也比较欠缺；同时很难使用复杂的速度函数描述界面的演化：根据单元密度 值抽取的界面位置也不十分精确。 v o f 方法将研究区域离散成网格单元，l e v e ls e t 方法则将区域离散成网格 点。l e v e l s e t 方法1 i 仅兼具v o f 方法的优点，而且还有许多优良的性质。l e v e ls e t 模型将运动界面表示为高一维场函数的零等值面，通过l e v e ls e t 方程驱动场函数 演化，从而跟踪界面演化的动态过程。l e v e ls e t 方程为一阶对流双曲方程，满足 双 l j 守恒律，通过缩小网格步长、应用高阶数值格式和使用自适应方法等手段， 可以提高数值解的精度，通过方程的初始值可以确定计算结果的误差范围；对界 面的法向、曲率等内在几何量的估算也十分简单而精确；可以自然地应用复杂的 速度函数。这许多优点使得l e v e ls e t 方法成为跟踪界面演化的重要方法。但l e v e l s e t 方法也存在不足之处，动态l e v e ls e t 方程依赖于时间参数( t i m e d e p e n d e n t p d e ) ，数值求解过程十分耗时：界面隐式化为场函数的零等值面，演化过程不 够直观；对界面细节的捕捉不敏感，不能保持初始界面的特征信息；从场函数中 抽取的界面也不够精确。总的来说，l e v e ls e t 模型是全局的方法，为界面演化、 变分问题的研究提供了一个行之有效的解决方案，许多实际问题也都可以归结为 设计合理的速度函数代入l e v e ls e t 模型求解。近年来的许多工作正是基于l e v e l s e t 方法的特性展开的。 i 2l e v e ls e t 方法的研究现状 o s h e r 与s e t h i a n 提山了l e v e ls e t 模型，同时给出了严格的数学推导过程和 数值求解格式 s e t 8 5 ，s e t 8 7 ，o s 8 8 1 ，为后续研究工作奠定了坚实的理论基础。在 众多研究人员的共同努力下，l e v e ls e t 方法得到卅i 断发展和完善。下面我们简要 介绍l e v e ls e t 方法在图形图像领域中一些具有代表性的研究成果 b s 0 2 ， 浙江人学坝i 学位论文“v e ls e t 方程的求解1j f 皿用 第一章绪论 b f m + 0 4 ，s e t 9 6 b ，o f 0 1 ，o f 0 2 。 1 2 1 原始模型的改进 在高一维空间中求解动态l e v e ls e t 方程是个耗时的过程，如何提高计算效 率对l e v e ls e t 方法的实际应用具有重大的意义。当界面运动的速度函数恒正( 负) 时，l e v e ls e t 方程转化为与时间无关的静态方程，s e t h i a n 提出了f a s tm a r c h i n g 方法快速求解 s e t 9 6 a ，s e t 9 6 b ：根据运动界面隐式化为场函数的零等值面这一性 质，窄带法将求解区域限定在界面所处区域周围的窄带内，减少了大量计算 【a s 9 5 ；s u s s m a n 等在求解二相流问题时提出自适应投影方法求解n a v i e r - s t o k e 方程，同时提出了自适应l e v e ls e t 方法提高捕捉流体界面的精度 s f s 0 9 8 ， s a b + 9 9 ：b a r t h 等提出了三角域上的l e v e ls e t 方法 b s 9 8 】，x u 等进一步提出了 三角域上的自适应方法，应用于图像分割【x 硎】。本文第三章将对各种自适应 方法做一个总结和分析，并针对三角网格曲面演化的特性给出我们的自适应求解 方法。l e v e ls e t 方程是h a m i l t o n j a c o b i 方程的一种具体形式，满足双曲守恒律， 研究人员还提出了许多高精度的数值计算格式，最具有代表性的有e n o 和 w e n o 格式等 o s 9 1 ，l o c 9 4 ，j p o o ；针对不同的研究领域和研究对象，也有许 多巧妙的改进方法。 1 2 2 网格模型的处理 计算机图形学通常使用三角网格表示物体。通过三维扫描仪等手段可以获取 实体模型采样点的空间位置信息，但这些原始数据必须通过重建、去噪、编辑等 处理后才能真正投入使用。l e v e ls e t 方法在网格模型的数字几何处理中起到了重 要作用。 通常散乱点网格重建是对模型进行数字几何处理的第一步。z h a o 等将重建的 i 角网格曲面定义为散乱点张成的极小 # 1 面，提出了变分l e v e ls e t 方法 z o m k 0 0 ， z o f 0 1 】。该方法首先在散乱点周围区域内构建初始网格，定义变分能量方程驱 动界面演化，利用l e v e l s e t 方法能有效处理拓扑结构改变的优点，得到了很好的 三角网格曲面重建结果，如图1 1 所示。 浙江人学坝1 。学位论文k v e ls e t 方程的求解i j 脚用 第一章绪论 l i n k e d t o r i 模型曲面重建 z o m k 0 0 b u d d h a 模型曲面重建 z o f 0 1 幽1 1 散乱点网格重建 t a s d i z e n 等将图像与信号处理的技术应用于网格处理。算法首先建赢叟分方 程修正顶点法向朝向，再根据新的法向信息改变顶点的空间位景，顺次迭代，得 到曲面去噪、光顺或增强的效果 t w b 0 0 2 】，如图1 2 所示。 图1 2 网格曲面去噪、光顺与增强 t w b 0 0 2 】 图1 3 使用l e v e ls e t 算子处理模_ _ j _ m b w b 0 2 】 4 浙江人学硕i 学位论文k v e ls e t 方程的求斛，j 城用第一审绪论 m u s e t h 等通过设计不同的速度函数定义出一系列的曲面编辑算子，对几何模 型进行去噪、编辑、雕刻等操作 m b w b 0 2 ，m b w + 0 5 ，如图1 3 所示。 从源物体到目标物体的m o r p h i n g 动丽是非常重要的特效。b r e e n 等使l e v e l s e t 方法生成不同亏格网格之间的m o r p h i n g 动画，中间曲面的拓扑结构变换过程 十分自然，如图1 4 n 示。他们还提出了稀疏场方法，提高方程求解的效率 【b m w m 0 1 ，b w 0 1 。 图1 4l e v e ls e t 方法应用 - m o r p h i n g 动画 b m w m 0 1 l e v e ls e t 方法还被应用于等距线生成 k b 9 3 ，极小曲面计算【c h 0 9 3 】，跟踪 晶体生长过程 s s 9 2 ，g f c 0 0 3 等。 1 2 3 图像视频处理 在图像视频处理中，经常需要进行图像增强、分割、去噪、修复等操作，特 别是在医学图像处理中，特征识别非常重要。l e v e l s e t 方法同样提供了强有力的 处理手段。二维空间中的l e v e ls e t 演化表现为活动轮廓线( a c t i v e c o n t o u r s ) 。 k a s s 等提出的s n a k e 方法是应用于图像分割的经典算法 k w t 8 7 ，c a s e l l e s 等 c k s s 9 7 和m a l l a d i 等 m s v 9 4 1 将l e v e ls e t 方法应用于图像分割，同样达到了 非常好的效果。算法首先在图像上任意给出初始曲线，通过每个象素的能量值、 梯度值等参数设计速度函数，驱动活动轮廓线演化，最终到达能量最小的位置， 如图1 5 所示。在此之后，研究人员还提出了许多更好的能量函数定义方法改进 卜述工作 s h a 9 6 ，c k s 9 7 ，p d 9 8 】，c h a n 等基于m u m f o r d s h a h 分割函数求解最小 分割问题的方法也非常经典 c v 0 1 】；图像处理中的许多方法还被应用于体数据可 见化 l k h w 0 4 、体数据分割 w b m s 0 1 ，o p 0 3 ，b w m z 0 5 幕f l 网格分割 y g z s 0 5 等领域。 浙江人学蜘 i 。学位论文l e v e ls e t 方程的求样j 随用鹕一章绪论 幽1 5l e v e ls e t 方法应用丁_ 图像分害i c k s s 9 7 1 b e r t a l m i o 将图像视为不可压缩流，通过求解n a v i e r _ s t o k e s 方程将信息传递致 损坏的区域进行图像修复( i n p a i n t i n g ) b b s 0 1 ；b a l l e s t e r 等求解图像关于灰度值和 梯度方向的变分差分方程，在信息丢失的区域光滑捅值，修复图像 b b c + 0 1 】， 如图1 6 所示。 图1 6 图像修复 b b c + 0 1 根据视频图像序列具有连续性的特性，s h i 等提出的实t c l e v e ls e t 方法可以 达到实时跟踪视频中物体运动边界线的效果 s k 0 5 ，如图1 7 所示。该方法使用 象素单元之间的交换操作来演化边界线，同时引入高斯滤波函数保持边界线的光 滑性，替代了费时的l e v e l s e t 方程的求解过程。 幽1 7 实时l e v e ls e t 方法跟踪视频中物体的运动轨迹 s k 0 5 】 1 2 4 流体动画 跟踪流体自由表面的演化是l e v e ls e t 模型最基本的应用领域。f e d k i w 等提出 的粒t l e v e ls e t 方法可以达到更高的精度，在表面上自适应分布的粒子不但被用 于修正距离场函数值，而且可以模拟水滴飞溅的效果。该方法已经成功应用于电 浙江人学删 j ：学位论文_ l e v e ls e t 方程的求解。，脚用 第一章绪论 影特效制作 f f 0 1 ，e m f 0 2 ，e f f m 0 2 ，e l f 0 5 】。近几年出现的许多基于物理的方法 模拟流体、火焰等动画效果，也都使用粒子l e v e ls e t 方法捕捉动态演化的界面 【f s j 0 1 ，c m b + 0 2 ，n f j 0 2 ，s t a 0 3 ，g b b 0 4 ，l g f 0 4 ，z b 0 5 ，h k 0 5 ，g s l f 0 5 ，如图1 8 所示。 ( a ) 自由表面 e m f 0 2 】( b ) 粘性流体 g b b 0 4 】( c ) 曲面上的流 s t a 0 3 】 图1 8l e v e ls e t 方法应h = i 于流体动画 1 3 本文工作 本文以l e v e ls e t 方程的求解与应用为中心，讨论了l e v e ls e t 方法的研究与 应用现状，针对一些具体问题和现有方法的不足之处，提出了新的解决方案。最 后结合作者的工作，对该领域今后的研究工作做了一个展望。 本文的主要贡献有以下几点： 提出了基于窄带的自适应l e v e ls e t 方法 针对三维空间中的三角网格曲面模型演化的特点，本文提出了基于窄带的自 适应l e v e ls e t 方法。自适应方法构建粗网格满足界面演化的整体需求，同时估 算粗网格点的曲率值，使用快速扩散法聚类高曲率点，估算聚类点集的最小包围 盒，以此为基础构建细网格捕捉界面演化的细节。粗、细网格均为独方的计算单 元，场函数在自适应网格中演化，从而跟踪界面演化的动态过程。我们给m 了一 套简单实用的数据结构，避免了多个网格之间频繁的坐标变换和捅值操作。最后 我们将白适应方法应用于三角网格模型的m o r p h i n g 动画。实验结果与误差分析 表明，自适应方法能有效减少计算量，达到更好的界面跟踪效果。 提出了混合等距线算法 本文结合网格方法与矢量方法的优点，提出了混合等距线算法。算法求解静 浙江人学坝f ：学位论文b v e 】s e t 方程的求解l j 成用 第一章绪论 态l e v e ls e t 方程生成距离场，从中抽取等距线的拓扑结构信息，通过顶点调整、 奇异处理、特征识别等步骤最终将等距线表示为基本图元的组合。我们以上述算 法为基础实现了一个原型系统，同时给出了系统流程和相应的数据结构。基于流 形上的静态l e v e ls e t 方程求解，本文还实现了带边界三角网格哺面上的广义等 距线和测地线算法，进而生成参数网格计算原始曲面的平面参数坐标，达到了较 好的参数化结果。 本文各个章节的内容安排如下： _ 第一章讨论l e v e ls e t 方法在各个主要领域中的研究与应用现状及其优点 和存在的问题，对当前的研究热点进行了归纳总结。 _ 第二章简要介绍l e v e ls e t 模型的理论背景，同时给出了经典的求解算法 和一些相关工作。 _ 第三章主要讨论动态l e v e ls e t 方程的求解算法，对现有的自适应方法进 行分析和总结，进而介绍基于窄带的自适应l e v e ls e t 方法。 一第四章讨论静态l e v e ls e t 方程的应用。主要介绍了混合等距线算法和带 边界三角网格的平面参数化算法。 _ 最后一章是对本文的总结，并提出了展望。 浙江人学坝l 。学位论文l e v ds e t 方程的求解4j 崩用 筘一章l e v e ls e t 方泣 2 1 引言 第二章l e v e ls e t 方法 l e v e ls e t 模型的推导与求解是进行深入研究的基础，本章介绍l e v e ls e t 模 型的理论背景和求解算法。 2 2 l e v e ls e t 模型 考虑维空间吼”中的一张m l 维光滑闭合超曲面r ，r 界定了吼“中的一个 开区域q ，f 上的点沿着该点法向以给定的速度函数，运动，f 可以是多个变量 的函数，如超曲面的曲率、法向、空间位置等，同时忽略f l f | 面的切向速度，如图 2 1 所示。用r ( f ) 表示r 运动得到的曲面簇，t 为时间变量。给定f 之后，目标就 是要构造r 的嵌入场函数，通过运动方程驱动场函数演化从而跟踪r 的演化 o s 8 8 ， s e i 9 6 b ，0 f 0 2 1 。 蚓2 1 界面r 以速度f 沿法向方向运动 s e t 9 6 b 若f 恒正( 负) ，则界面将总是向外( 内) 扩张，r 的运动小依赖j - , 时间变 量t ，到达空间定点的时间函数值r ( x ，y ) 唯一确定，如图2 2 所示。 浙江人学倾i j 学位论文l e v e ls e t 方程的求解i i 应用 第二章l e v e ls e t 方浊 f0 、 、 -f 7， l ， ，一 ) ( f1 ， t ( x - 图2 2 曲线以述度f 0 演化剑达点“y ) 的时间兀s e t 9 6 b 】 位移与速度的关系为l = f d ，i d 。其中d t 为时间微分，出为d t 时间内的位 移，界面上的点在f 的作用下沿时间函数r 的梯度方向v7 1 运动，则 f i v t ( x ) if = l i r = f r ( x ) ，= o ( 2 1 ) 上述方程即为静态l e v e ls e t 方程。如果速度函数f 的设定仅依赖于r 的空间 位置，则方程( 2 1 ) 即为e i k o n a l 方程 s e i 9 6 a ，s e i 9 6 b 】。 如果f 不是恒正( 负) 函数，则界面到达空间定点的时间函数t ( x ，y ) 是个 多值函数，此时构造新的嵌入场函数妒，硇满足： f ( t ) = xi 妒( x ，t ) = 0 ) ( 2 2 ) 其中x ( f ) 表示点x 随时间变化的轨迹函数，x ( r = o ) 即为初始曲面l ( w - 0 ) 上的 点。场函数旷般取为有向距离函数 条件 x e q x a q = f ( t 1 x 芒q ( 2 3 ) 其巾d ( x ) = m i n ( i ip x1 1 ) 表示点x 到曲面r ( f ) 自日最近欧氏距离值，驴满足边界 v p e f ( f 1 f ( t = 0 ) = x1 0 ( x ，t = 0 ) = 0 ) ( 2 4 ) 由于点x 沿着法向i 方向运动，得x ，月= f ( x ( f ) ) 。方程妒( x ( f ) ，f ) = o 两边 分别对f 求导： i o 掣 r，【 = ) x 浙江人学帧 。学位论义k v e ls e ! 方程的求解j 应用 第= 市l e v e ls e l 方浊 谚十v 妒( x ( f ) ，t ) x ( f ) = 0( 2 5 ) 再由x ，亓= ，( x ( r ) ) ，法向f i = v 妒i v 妒i 和边界条件可得： 痧,钏+fivof糍0 l ( 2 。) i = 矿( j ，) i f _ l 、。 上述方程即为动态l e v e ls e t 方程。 无论是静态方程或动态方程，相应的嵌入场函数均为单值函数，运动界面r ( f ) 在任意时刻的空间位置由t ( x ，y ) = r 或者妒( x ，y ，r ) = 0 给出。若场函数以连续的速 度函数演化，则r 始终保持为场函数的零等值面。在高一维空间中演化场函数， 可以鲁棒而自然地处埋拓扑分裂或合并、自相交等复杂的情况，如图2 3 所示。 。，1h rl e v e l t 删恕仲母( d a r k a ，1 ： h l r ；f lt i m e ：i 抽h d e d 暾印。慷yh a s 柏a n s t m l ， 1 、州) 蝌l 瑚m 抽i n i l i a ln l 慷( i nl i 曲g r a y ) 一 y i d d i n g8s i f l g l c u r v e * i nt h e ；赫l e v e l 蝌t 圈2 3 嵌入场函数在各个时间片的零等值厩 s e t 9 6 b 2 3l e v e ls e t 方程求解 静态与动态l e v e ls e t 方程可以统一写为h a m i l t o n j a c o b i 方程： 伽，十h ( d u ，x ) = 0( 2 7 ) 其中d u 表示对“各个分量的偏导，如，u y ，u ：。动态l e v e ls e t 方程取a = l ； 静态l e v e ls e t 方程取a = 0 ，h = f i v “i 一1 。方程( 2 7 ) 满足双曲守恒律，可以通 过相应的逆风差分格式求得满足粘性条件的弱解。 2 3 1 静态方程求解 。 静态l e v e ls e t 方程的特征线为币直于界面向外的直线，方程的解总是沿着 特征线方向传播。利用这一特性，s e t h i a n 提出了f a s tm a r c h i n g 算法。该算法使 浙江人学坝卜学位论文k v ds e t 方程的求斛i ，应用帮二章l e v e ls e t 方泣 用最小堆( m i nh e a p ) 数据结构，每次弹出当前活动点集中具有最小函数值的网格 点，求解相邻网格点的函数值 s e t 9 6 a 。在孵中算法的具体步骤如下： 1 初始化 a ) 初始网格点标记为活动点( a l i v e ) b ) 其余网格点标记为无穷远点( f a r a w a y ) ，函数值设为无穷大 2 扩散： a ) 开始循环：从活动点集中取出具有最小函数值的网格点p ( x l ，y ，z 。) ， 标记为固定点( s t a b l e ) b ) 分别考虑p 点的邻近网格点p ( 十y j , z ) ， 尸( ”y ，z 。) ， p ( ，yj _ l ，z k ) ，p ( x j ，y j + i z t ) ，p ( ，y j ，z k - i ) ，p ( ，y ，z k + 1 ) ： i 若为固定点，则继续，否则 讧 标记为活动点，根据差分格式计算该点函数值，比较当前函数值， 取绝对值较小者 c ) 开始新的循环，直到活动点集为空时退出 f a s tm a r c h i n g 算法的执行过程如图2 4 年1 1 2 5 所示。根据速度函数的设定，可 以直接计算界面所处区域周围的网格点值作为方程求解的初始值。若曲面以单位 速度演化，则到达时间函数t ( x ，y ) 与距离函数值相等。 o0 000 o( ro0o0000g 。 霉 ：。 ” ； 啦。 + ” ) 。一- 。 。_ _ _ l - 卜 。 i 一日- _ 卜 。 - i； 。7。、+。：3。；1啦e 。+ 一+ t3 ， 、，：j 。 ：， ? ，+ 。 ，)! ，j 、 图2 4f a s tm a r c h i n g 算法执行过程 s e t 9 6 b 幽2 5f a s tm a r c h i n g 过程并个网格点的状态 s e t 9 6 b 】 浙江人学坝1 峙位论史l e v e ls e t 方程的求解，应用 第一章l e v e ls e t 方法 静态l e v e ls e t 方程( 2 1 ) 在吼3 中的满足粘性条件的一阶逆风差分格式为： 陋m i n ( 啄d i i t t 7 浆+ n 心m 眩( d s r t 黧删r a i n ( 端巍：偿s ， ，o ) ： ，o ) 。+ d 荽r ，o ) 。】= l ，2 由于时间函数的特征线垂直于界面，上述格式等价于： m a x ( m a x ( q - ，。* t ，o ) ，一n f m ( o ；t t ，o ) ) 2 + m a x ( m a x ( p 嘉t ，o ) ，一m i n ( d i 2 t ，o ) ) 2 + ( 2 9 ) m a x ( m a x ( d i 7 1 ，o ) ，一m i n ( d 茹r ，o ) ) 2 = 1 ! 不妨设到达时间函数值为正，f o ，同时有 i = m a x ( m 呱d 蔷r ，o ) ，- m i u ( d i i , t ，o ) ) 2 7 2 = m a x ( m a x ( d i ； t ，o ) ，一i i l i n ( 口丁，o ) ) 2 五= m a x ( m a x ( d ；i = 7 _ ，o ) ，一i t l i n ( d 茹7 1 ，o ) ) 2 满足n 疋乃 其中喀，d y ，4 分别为z ，y ，z 方向的网格步长，令a = l d , 2 ，b = l d ，2 c ：1 d 2 ，未知量，为所求网格点的时间函数值。 由于f a s t m a r c h i n g 方法在每一个计算步骤总是弹出当前具有最小值孔的网 格点，且距离越远的网格点，到达时间函数值越大，因此有：t i t 。 若：t l r 死死，则为单点情形： t = 五十以 ( 2 1 0 ) 否则，若t 1t 2 t t 3 ，则为两点情形： r ：至兰圣竺地竺) 二塑圣二壁 ( 2 1 1 ) a + 占 否则，若t j 疋乃 7 1 ，则为三点情形： r ：互生型堑巫互亟亟豆匝亟噩噩画( 2 ，1 2 ) a + b + c 通过公式( 1 0 ) 一( 2 1 2 ) 计算得到的t 值必然大于7 i ，也保证了f a s tm a r c h i n g 算法的正确性和有效性。 浙江人学颂i + 学位论立l e v e ls e t 方程的求斛j 舸用 第一章l e v e ls e t 方法 2 3 2 动态方程求解 给定凸速度函数f ，孵中的动态l e v e ls e t 方程( 2 6 ) 可以用下列一阶逆风差 分格式求得满足粘性条件的数值解： 蠕1 = 谚n m a t e m a x ( f , i k ，o ) v + + m i n ( f j o ) v 一 ( 2 1 3 ) 其中 v + = m 一。- ao ) 2 + m i w , 7 ，o ) 2 + m a x ( d j ，o ) 2 + m i n ( d j ，o ) 2 + m a x ( d 叠，o ) 2 + m i n ( d 蔷，o ) 2 】，2 v 一： n l i n ( 聪 0 ) ：+ m a x ( 菇，o ) ：+ l t l i n ( 蠢，0 ) 2 + 但1 4 m a x ( p ，o ) 2 + m i n ( d i f ，o ) 2 + m a ) 【( 磙，o ) 2 r d 茹( 啄) 为聪妒( d i 妒) 的缩写，其余记号的含义类似。 高阶格式及速度函数非凸的情况比较复杂，本文暂不作论述，详细阐述见文 献 p e n 9 0 0 s e i 9 6 b 。 曲面法向与中曲率的估算公式为： 法向： 拈尚_ c 而寿酽，而者酽，而孝酽 中曲率为单位法向的散度： f ( 嘞+ 缱：) 破2 + ( 破：+ 如) 砖2 + ( 屯+ 略) 晓2 1 瑚尚= 也塑器蒜掣 上述公式对静态方程同样适用。 在整个网格中求解动态l e v e ls e t 方程的时间复杂度为o ( n 3 ) ，效率十分低 下。由于界面隐式化为场函数的零等值面，一般只需要求解场函数在界面周围区 域中的网格点值，距离较远的网格点对计算结果的影响可以忽略不计，窄带法将 力。程的求解区域限定在界面周围的一斟窄带内，时间复杂度降为o ( 1 c 9 2 ) ，节省 了大量的计算 a s 9 5 。算法的具体步骤如卜- ： 1 设定窄带包含的网格层数阈值 2 窄带内部的网格点标记为活动点( a l i v e ) ：边界上的网格点标记为边界点 1 4 浙江人学硕i 。学位论义h v e ls e t 方程的求斛l i 一用 第一章l e v e l s e t 方法 ( l a n d m i n e ) 其余网格点标记为无穷远点( f a r a w a y ) 。若网格点位于界面 外部( 内部) 则函数值设为无穷大( 小) 3 设定初始值，按动态方程求解格式计算窄带中网格点的函数值。若界面 到达窄带的边缘区域，则构建新的窄带 4 开始新的循环 招主过 一中。中一+ 。 l ，女- li 寸_ t 捌 留丢? x”广? 捌 军h 豫翮 忖t q r ”、 r 倒2 6 窄带法演化界面 s e t 9 6 b 使用窄带演化界面的过程如图2 6 所示，求解动态方程的初始值可以通过将 速度函数取为单位值，使用f a s t m a r c h i n g 方法生成。时间步长血由该时刻窄带内 网格点上的最大速度值决定，满足c f l 条件：a t 0 xeq 口 ( 3 5 ) ( 工) 0 工匹q b b a 图3 ，1 5 源物体演化至目标物体 b m w m 0 1 当两张曲面部分重叠时，驱动q 。内部的源曲面膨胀成s 。，同时q 。外部 的源曲面收缩到上，当s = s 。时达到平衡，演化过程结束，如图3 1 5 所示。 满足公式3 5 的驱动函数虼一般通过为两个曲面间有向距离场函数的差值计算 得到。我们使用基于窄带的自适应l e v e ls e t 方法，得到了自适应光滑变形的效 果。 定义源曲面为一平面，目标曲面为m a x p l a n c k 曲面模型。自适应方法在演化 过程中自动检测拂面的曲率值，构建自适应细网格捕捉逐渐出现的细节，提高汁 算精度。图3 1 6 展示r 源曲面演化成为目标曲面的过程和构建的自适应绌网格及 旧格点。 浙汀人学蛳i j 学位论文k v ds e t 方程的求解j 应用第二三章罐于窄带的白适应l e v e ls e t 方法 ( a ) 迭代5 0 个时间步( b ) 迭代4 0 0 个时间步( c ) 迭代7 5 0 个时间步 图3 1 6 白适应l e v e ls e t