基于近似几何误差的动态隐式曲线重构.doc

1、Vol.15, No.62004 Journal of Software 软 件 学 报1000-9825/2004/15(06)0960基于近似几何误差的动态隐式曲线重构 杨 周旺 +, 邓 建松 , 陈发 来(中国科学技术大学 数学系, 安徽合肥 230026)Dynamic Implicit Curve Reconstruction Based on Approximate Geometric Distance*YANG Zhou-Wang+, DENG Jian-Song, CHEN Fa-Lai(Department of Mathematics, University of Sci

2、ence and Technology of China, Hefei 230026, China)+ Corresponding author: Phn: +86-551-3601009, Fax: +86-551-3601005, E-mail: , http:/Received 2004-06-10; Accepted 2004-07-05Yang ZW, Deng JS, Chen FL. Dynamic implicit curve reconstruction based on approximate geometric distance. Journal of Software,

3、 2004,15(6):960974.http:/ An implicit curve reconstruction method is proposed which represents the curve with an algebraic tensor-product B-spline, and minimizes the tension of the B-spline and the approximate geometric distance between the curve and the point set. The method is dynamic and self-ada

4、ptive based on trust-region algorithm in optimization theory. The specification of the initial shape is priceless, and the high-quality reconstruction curve is obtained in a robust way. Some examples are implemented.Key words: Implicit reconstruction; geometric error; algebraic tensor-product B-spli

5、ne; trust region; dynamic摘 要: 本文提出一种以代数张量积样条曲线作为几何表示形式，基于近似几何 误差和薄板能量极小化的隐式重构模型。同时结 合最优化理论中的信 赖域思想， 给出自适应的迭代求解算法及其 实现。这种方法采取无代价初始化技术，通过迭代能 稳定地达到目 标点集的高质量重构，特 别是对 复杂形状的目标，具有很强的处理能力。关键词: 隐式重构;几何误差;代数张量积样条曲线;信赖域;动态方法中图法分类号: TP391 文献标识码: A1 引言在 CAD/CAM、计算机图形学、科学计算可视化和图像处理等领域以及实际工业部门中经常碰到这样的问 Supported b

6、y the National Science Fund for Distinguished Young Scholars under Grant No. 60225002 (国 家 杰 出 青 年 基 金 ); the Teaching and Research Award Program for Outstanding Young Teachers in Higher Education Institutions of MOE (教 育 部 高 校 优秀青年教师教学科研奖励计划); the National Research Foundation for the Doctoral Progr

7、am of Higher Education of China under Grant No. 20010358003 (教育部博士点基金); the National Natural Science Foundation of China under Grant No.10201030 (国家自然科学基金 ); the Research Foundation for Young Scholars in University of Science and Technology of China under Grant No.kb0122 (中国科学技术大学青年基金)作者简介: 杨周旺(1974

8、),男,福建福安人,讲师,主要研究领域为隐式曲线曲面造型; 邓建松(1971 ),男,博士,副教授,主要研究领域为计算机辅助几何设计与计算机图形学;陈发来 (1966),男,教授,博士生导师,主要研究领域为计算机辅助几何设计与计算机图形学.杨周旺 等:基于近似几何误差的动态隐式曲线重构 961题：对给定的或采集到的表示某一几何形状的无序点集（又称点云） ，如何重构曲线（面） ，使得目标点集离曲线（面）在某种度量下偏差最小。这就是经典的曲线和曲面拟合问题。但是随着计算机技术和数据采集技术的发展，所得数据量非常庞大，目标曲线（面）的拓扑结构相当复杂，因此对这类问题的研究显得十分困难，近年来已形成了

9、一个新兴的研究领域。这种重构技术为产品的快速开发和原型化设计提供了有效的途径，它可应用于机械、轻工、汽车、航空以及科学计算可视化、医学图象处理等领域。在曲线（面）重构中，可以根据所得几何形状的表示形式，把重构方法分为参数曲线（面）重构和隐式曲线（面）重构。另外，由于所给点云通常具有噪音（noise） ，因此最后重构出来的曲线（面）一般是拟合给定数据点。为了建立起稳定的重构算法，通常采用迭代的方法，即从某个初始形状出发，通过逐步减少某种度量，使得最终形状收敛到度量极小的状态。这就是动态重构方法。本文给出的算法是一种基于隐式表示的动态重构方法。我们采用代数张量积样条曲线表示结果。从上世纪八十年代开

10、始，已有大量研究者在参数曲线和曲面重构方面进行了研究。例如，Kass, Witkin 和Terzopoulos1用参数型的动态轮廓(active contours, 亦称 snakes)拟合图像的边界； Blake 和 Isard2使用参数样条曲线技术，以法向量偏差为误差项作动态逼近；Pottmann 等 3,4研究了动曲线上的采样点到目标曲线的平方距离，通过求解一个优化问题，动态调整样条曲线；Wang 等人 5,6改进了 Pottmann 的结果，并对方法进行了理论分析工作。而在隐式曲线和曲面重构方面，Pratt 7给出在规范约束下代数误差极小化的拟合方法，以及 Taubin8的带权代数误差

11、极小化模型，用于目标物体的识别与表面分割。由 Osher 和 Sethian911提出的水平集(level sets)方法，用于解决临界面的追踪问题，这是一种隐式的动态轮廓方法；Zhao 12通过建立能量模型，应用基于变分的水平集方法进行重构；Foster 和 Fedkiw13综合物理模型与水平集方法建立水的模型等等。在上述隐式曲线和曲面重构算法中，得到的结果是离散表示。 Jttler14,15先估计目标数据点集的法向信息，综合考虑代数误差和法向偏差，建立基于代数张量积样条函数表示的隐式重构模型，最终通过求解线性方程组得到重构曲线或曲面。在用参数曲线（面）进行重构时，需要对无序数据点进行参数化

12、，这并不是一个简单的问题。而且利用参数表达难以构造出具有复杂拓扑结构的曲线（面） 。另外，在动态参数曲线（面）重构中需要初始形状与所构造的曲线（面）具有大致相同的拓扑结构。因此参数曲线（面）重构方法往往很难满足基于点云数据的曲线（面）重构的设计要求。与此相对，采用隐式曲线（面）表示进行重构时，无需对点集进行参数化，而且可以构造出具有复杂拓扑结构的曲线（面） 。通常在动态隐式曲线（面）重构中，迭代过程对于初值的敏感性并不是很高。实际上，参数曲线（面）和隐式曲线（面）(通常采用的是代数曲线（面） ，即多项式的零点集) 是计算机辅助几何设计中表示曲线（面）的两种主要方式。隐式曲线（面）除了具有上述相

13、对于参数曲线（面）的特点外，它还具有其它的一些优势。例如，可以很容易判断点是否在隐式曲线（面）内外。关于参数曲线（面）和隐式曲线（面）的详细讨论，请见 16,17。本文着重在于构造一个有效的隐式曲线重构算法，因此我们只讨论从平面点云构造曲线的问题。所给的方法很容易推广到从三维点云构造曲面的情形。本文的组织如下：第二节描述点云到隐式曲线的几何误差度量以及近似表示；第三节介绍一种特殊的隐式曲线，即代数张量积样条曲线，及其几何特征；第四节提出近似几何误差和薄板能量极小化的隐式重构模型并给出求解的自适应迭代算法；第五节则是具体实现和一些重构例子；第六节概述了本文的主要结果并展望进一步的工作。2 几何误

14、差在本节我们首先阐述了隐式曲线重构问题，然后给出了从隐式曲线到点云的一种几何误差度量，以及这种度量的近似表示。2.1 隐式曲线重构问题设 是区域 上的光滑函数（具有一阶连续偏导） ，那么可由隐式方程 定义 的零RfK: 0)(Xff962 Journal of Software 软件学报 2004,15(6) 点集(2.1)0)(|)( XfRXfVK当 时 一般是平面曲线，当 时 一般是曲面。2K)(fV3K)(f本文我们只讨论 时的情形。给定二元函数空间 以及平面上的一组点 ，用2KSNiX1表示零点集 到 的一种偏差。那么隐式曲线重构问题就表述为)(,(Er1fXNi )(fVNiX1(

15、2.2)(,(Er1mnfViSf 对于具体的 ，我们称 为 到 的误差。Sf),(Er1NifNiX12.2 误差的定义定义 2.1. 固定 以及一点 ，用 表示从 到 的欧氏距离。定义fX)(,fVd)(fV(2.3)NiiNig fXdf121,Er为几何距离误差。其中 表示 到 的几何距离。YfVdfVY)(mn)(, )(f由于 的非线性， 一般没有解析表达，因此利用几何距离误差进行隐式曲线重构是非常难以求解)(fVX的。Sampson 18在讨论二次曲线的拟合时，提出几何距离误差的一阶近似，(2.4)NiiNis XffVX121)()(,(Er以下我们将这种近似几何误差称作 Sa

16、mpson 误差。另外，也可以采用如下方式定义误差：(2.5)NiiNiaffVX121)()(,(Er我们称之为代数误差。Bajaj 等 19曾基于这种误差讨论了应用隐式代数曲面进行高阶插值与拟合的问题。与前面两种误差的定义相比，这种定义不具有坐标变换不变性。但如果基于它进行曲线重构，对应的问题是一个线性最小二乘问题。几何距离 Sampson 距离 代数距离图 1 三种距离的比较 图 2 矩形框为点云法向不存在的区域杨周旺 等:基于近似几何误差的动态隐式曲线重构 963图 1 给出了在三种距离意义下到代数曲线 的等距曲线图。0xy3 代数张量积样条曲线双 次张量积样条函数的定义为d(3.1)

17、srsryNxMpyxfXf, )(),()其中 是控制系数。基函数 分别是对应于节点序列 和.,1,nsmrRprs )(,sr 1dmr的 次 B-样条 20。函数 的零点集 称为代数张量积样条曲线。1dns)(Xf)(fV在隐式曲线重构中，取适当的矩形区域 为包含目标数据点集 并略微扩大,11dndm NiX1的区域，且令 是满足 的内部等距节点序列。,11d , , 1dnn本文将应用这种曲线进行动态隐式重构。3.1 代数张量积样条曲线的几何特征基于上面表示的代数张量积样条曲线，我们可以很方便地计算给定目标点集 上的代数距离和梯度信NiX1息，以及定义在 区域上的薄板能量。记 为张量积

18、样条函数 的控制系数拉直向量，其中 ， ，),(1Lc),(yxf mnLrsnrjp)1(c。那么数据点 到隐式曲线 的代数距离定义为j, iX0(3.2)cisrisiri qyNxMpXf , )()(它是关于控制系数的线性函数。样条函数 在数据点 处的梯度向量是yxfiX(3.3)ciisirsi vuyNxpXf )()(于是，可用控制系数向量 来表示 Sampson 误差c(3.4)NiisBA1)(Ercc其中 都是 阶对称半正定阵。iiiii vuBqA, L张量积样条函数 在区域 的薄板能量定义为),(yxf(3.5)ccHdxyffTxy)2()2其中 （显然是半正定阵）的

19、元素对应于样条基 及其相应阶导数在区域 上的积分。根据样条H (),NMsr 964 Journal of Software 软件学报 2004,15(6) 性质，能解析求得 ，也可通过数值积分近似求得。H4 隐式重构模型与算法近年来在隐式曲线曲面重构方面，基于水平集方法如 12,13，通过解 PDE 得到离散表示的结果，往往计算量巨大。基于代数误差和法向估计的重构方法如 14,15，很大程度上依赖于事先对点云的法向估计，而估计法向本身是困难的，因为在某些情况下法向是奇异的（如图 2） 。本节所提出的动态隐式重构模型中，不需要事先估计点云数据的法向信息，在法向奇异处其重构效果也很好。4.1 优

20、化目标基于代数张量积样条曲线的表示，易得相应 Sampson 误差极小化的重构模型，但在实际问题中往往还应考虑重构曲线的光滑性。而薄板能量 反映了曲线的硬度， 并在一定程度上影响重构曲线的光滑性。在这里我)(cT们给出如下隐式重构模型：(4.1)(Er21)(mincwTRs其中 是光滑项权重系数，一般取很小的正数。w由于重构模型(4.1)中的优化目标函数 带有分式，其最优解无法直接求出。我们将结合最优化理论中的)(c信赖域 21思想，给出一种自适应的迭代求解算法。4.2 信赖域方法信赖域法是最优化理论中求解非线性优化问题的一类重要方法，其基本思想是在迭代点 附近用一个二)(kc次函数逼近原优

21、化问题（通常称作逼近子问题）并构造试探步 ，且要求试探步在信赖域之内，即每次迭代)(kd时有一正数 （信赖域半径）使得 ，其中 是某一范数。那么，对于逼近子问题而言，信赖域法的keke)(d试探步 使得 是以 为中心的广义球 上最优的点。由于试探步具有这一性质，信)(d)()(kdc)(kc )(kkec赖域法不进行一维搜索，而是当试探步 满足某种下降条件时令 ，否则 ，然后再)(k )()()1(kdc )()1(kkc以某种机制给出下一步的信赖域半径 。1ke对于我们的隐式曲线重构模型(4.1)而言，从适当初始点 出发（令 ） ，考虑在迭代点 附近用二次)0(c0k)(kc函数(4.2)d

22、cddcd kkkk GgRwHJJwHJRQ 21)()(21)()() 4321)( 逼近 。其中(c Ni kiiikNikiNikiiNiki BAJBAJBAJBAJ 13)()(412)(312)()(21)( cccc 在这里，我们记(4.3)(21(kkwHJg(4.4)Gk43)杨周旺 等:基于近似几何误差的动态隐式曲线重构 965显然， 是目标函数 在 处的梯度向量， （半正定） ，它是对目标函数 Hessian 矩阵的一种修正。kg)(cR)(k 0kG为了使算法具有理想的总体收敛性，先确定一个步长上限 （即信赖域半径） ，并由此定义 的信赖邻域ke )(kc。假定在这个

23、邻域里二次模型是目标函数 的一个合适的模拟，然后利用二次模型)(kkec )(cR确定修正量 。事实上，每一步迭代等价于求解信赖域子问题：)(dQ)(d(4.5)2)()(s.tmin21k kkkeGgRQddc 由约束最优化的 Kuhn-Tucker 条件 22知，该子问题的解可由线性方程组 来表征， 即确定一个kkg)(，使得 正定，并求相应线性方程组的解 。0kkG )(kd另外，由于第 步实际下降量与预计下降量的比值：(4.6)()()( kkkk QRc很好地衡量了二次模型近似目标函数的程度。据此可以给出一个更简洁的模式，即通过自适应地改变参数（等效于选择信赖域半径 ） ，保持二次

24、模型与目标函数尽可能一致的同时加速收敛。如取适当的kke，当 时，令 ；当 时，令 ；其它情形令 。10211k k412k kk21 k14.3 算法框架结合前面的讨论，本小节给出一个隐式重构的算法框架，并将指出它的一些收敛性质。隐式重构算法步骤如下：初始步给定： .0:,10,21)( klowc第 步迭代格式为：k1) 给出 ，计算 ，若 停止。k,)( kGg,k2) 采用对称 QR 分解求出线性方程 的解 。kkgd)( )(k3) 由（4.6）式计算比值 。k4) 更新 参数：若 ，令 ；若 ，令 ；否则，令 。1kk412k ),21max(1lowkk k15) 若 ，置 ；否

25、则令 ，计算 ，并置 。0k)()(kkc )()(dcciBNc)1(6) 令 ，转 1)。1:由信赖域方法的收敛性分析 21知，上述隐式重构算法具有总体收敛性。即由算法生成的控制系数序列如果有界，则必存在聚点 满足一阶必要条件 。进一步，如果 在该聚点 处的 Hessian 矩)(kc )(c0g)(cR)(c阵是正定的（即成立最优性二阶充分条件） ，则 二次收敛于 。)(kc)(c966 Journal of Software 软件学报 2004,15(6) 在算法的迭代过程中，需要注意的是，样条函数 在数据点 处的梯度模 在分母中),(yxfiXciiBXf2)(出现。为了防止溢出，在

26、实际计算中当 时，取一正小数 （如 ）代替。0)(iXf 8105 具体实现与重构例子用代数张量积样条曲线实现目标点集的隐式重构，首先要读入目标点集数据，指定样条次数 ，并如第 3d节中所述给定包含所有点的适当矩形区域 及内部等距节点 ， 然后生成相应的张量111dndm积样条基函数 。这样就可以事先计算区域 上的薄板能量矩阵 （见(3.5)式） ， 以及目标点)(yNxMsr H处的对应信息向量 。iX,iivuq第 4 节中所给的重构模型求解算法，要求给定某一 （初始曲线）开始迭代。在这里，我们给出一种无)0(c代价初始化技术，即令(5.1).,1,)1()(0)0(1)( Ljpnsmr

27、prsnrj c其中 是某一常数，使得控制系数向量 确定一条包住所有目标点的样条曲线，并将它作为初始曲线。以上0p )0(是迭代求解之前的预处理过程。我们还认为，重构模型目标函数中薄板能量项的权重 在迭代求解过程中也应该是一个动态变化的量。即w当远离目标点集时取较大的权重来加强迭代曲线的硬度，保证曲线整体移动的一致性；而当靠近目标点集时取较小的权重使得迭代曲线在局部有更强的适应性，提高逼近效果。我们采取线性递减的简单模式来动态地改变权重，即设定一个合适的初始权重 ，每步迭代令0w(5.2).1,0(,:1kk在从图 3 到图 7 的重构实例中，采用双二次 B-样条，矩形区域 ，节点序列取为10

28、,（内部等距） ，权重调节比例 ，停机准则 （如取 ）10,8,10, 75.0kg210及 。图中的离散点表示给定的目标数据点集，最外层的虚线是根据前面介绍的无代价初始化技7.,3.21术给出的初始曲线，实线则是当前得到的迭代曲线。从实例一（图 3）和实例二（4）可以看出，这种隐式重构方法能够很好地刻画目标的尖锐特征；实例三（图 5）和实例四（图 6）有力地说明了方法对复杂拓扑结构曲线重构的有效性，而采用参数表达的重构方法在曲线拓扑结构未预知的情况下是难以做到的；实例五（图 7）则显示了对目标点集具有厚度的点云数据进行隐式重构的效果。在表 1 中列出了这些实例的计算信息，其中最后一列表示在

29、PIV-2.2G 微机上基于 Matlab 6.5 实现前述算法的运行时间。表 1 重构实例的计算信息重构曲线目标点数 kg迭代次数 sEr运行时间秒实例一 41 0.0018 30 0.0035 0.5实例二 187 0.0084 33 0.0071 2实例三 372 0.0100 39 0.0272 6杨周旺 等:基于近似几何误差的动态隐式曲线重构 967实例四 105 0.0059 49 0.1032 2实例五 969 0.0515 50 4.9043 256 结论文中提出的基于代数张量积样条曲线表示的隐式重构方法，从一无代价初始化曲线出发，通过极小化近似几何误差来不断修正迭代曲线，逐步

30、稳定地达到对目标点集的高质量重构。这是一种动态的隐式重构方法，它的一个突出优点是可以忽略点云的拓扑。该重构算法稳定，具有总体收敛性，对初值不敏感。第 5 节中实现的重构实例进一步验证了上述特点，同时反映出该方法收敛速度较快，重构效果好。本文给定的张量积 B-样条节点序列是等距的，另外，根据点云数据如何确定和生成具有自适应性的节点向量将是今后进一步需要分析的工作。这种平面曲线的隐式重构方法可以很自然地推广到三维空间曲面的情形。我们将在其它论文中对此进行讨论。a. 初始曲线和迭代曲线(k=5) b. 重构曲线(k=30)图 3 重构实例一a. 初始曲线和迭代曲线(k=5) b. 重构曲线(k=33

31、)图 4 重构实例二968 Journal of Software 软件学报 2004,15(6) a. 初始曲线和迭代曲线(k=2) b. 重构曲线(k=39)图 5 重构实例三a. 初始曲线和迭代曲线(k=3) b. 迭代曲线(k=6)c. 迭代曲线(k=18) d. 重构曲线(k=49)图 6 重构实例四杨周旺 等:基于近似几何误差的动态隐式曲线重构 969a. 初始曲线和迭代曲线(k=8) b. 重构曲线(k=50)图 7 重构实例五References:1 M. Kass, A. Witkin, D. Terzopoulos, Snakes: active contour model

32、s. International Journal of Computer Vision, 1988, 321-331.2 A. Blake, M. Isard, Active Contours. New York: Springer-Verlag, 1998.3 H. Pottmann, M. Hofer, Geometry of the squared distance function to curves and surfaces. Technical Report of Institute of Geometry, Vienna University of Technology, 200

33、2.4 H.Pottmann, S.Leopoldseder and M. Hofer, Approximation with active B-spline curves and surfaces. Proceeding of Pacific Graphics 2002, IEEE Press, 8-25.5 W. Wang, et al., Active B-spline curve: Using squared distance minimization for curve reconstruction from unorganized data points. Technical Re

34、port of Computer Science Department, Hong Kong University, 2003.6 H.P. Yang, W. Wang and J.G. Sun, Control point adjustment for B-spline curve approximation. Computer-Aided Design, in press.7 V. Pratt, Direct least-squares fitting of algebraic surface. ACM Computer Graphics, 1987, 21(4), 145-152.8 R

35、. Taubin, Estimation of planar curves, surfaces, and non-planar space curves defined by implicit equations with applications to edge and range image segmentation. IEEE Trans. Pattern Anal. Mach. Intelligence, 1991, 13, 1115-1138.9 S. Osher, J. Sethian, Fronts propagating with curvature dependent spe

36、ed, algorithms based on a Hamilton-Jocobi formulation. Journal of Computational Physics, 1988, 79, 12-49.10 J.A. Sethian, Level Set Methods and Fast Marching Methods. Cambridge University Press, 1999.11 S. Osher, R. Fedekiw, Level Set Methods and Dynamic Implicit Surfaces. New York: Springer-Verlag,

37、 2003.12 H.K. Zhao, S. Osher, B. Merriman and M. Kang, Implicit and nonparametric shape reconstruction from unorganized data using a variational level set method. Computer Vision and Image Understanding, 2000, 80, 295-314.13 N. Foster, R. Fedkiw, Practical animation of liquids. Proceedings of SIGGRA

38、PH01, 15-22.14 B. Jttler, A. Felis, Least-squares fitting of algebraic spline surfaces. Advances in Computational Mathematics, 2002, 17, 135-152.15 B. Jttler, Approximate implicitization via curve fitting. Eurographics Symposium on Geometry Processing, 2003.16 G. Farin, Curves and Surfaces for CAGD-

39、A Practical Guide, 5th Edition. Morgan Kaufmann Publishers, 2002.17 J. Bloomenthal, et al., Introduction to Implicit Surfaces. Morgan Kaufmann Publishers, Inc., 1998.18 P.D. Sampson, Fitting conic sections to very scattered data: An iterative refinement of the Bookstein algorithm. Computer Graphics

40、and Image Processing, 1982, 18, 97-108.19 C. Bajaj, I. Ihm, and J. Warren, Higher-order interpolation and least-squares approximation using implicit algebraic surfaces. ACM Transactions on Graphics, 1993, 12(4), 327-347.20 Carl De Boor, A Practical Guide to Splines. New York: Springer-Verlag, 2001.21 M.J.D. Powell, On the global convergence of trust region algorithms for unconstrained optimization. Mathematical Programming, 1984, 29, 297-303.