196 2014，50（13） Computer Engineering and Applications 计算机工程与应用 三角网格曲面上离散曲率改进算法 马献德 1，谢小峰 2 MA Xiande1, XIE Xiaofeng2 1.南京电子技术研究所 一部，南京 210039 2.南京禄口国际机场 动力技术部，南京 211113 1.Department of Overall Research, Nanjing Research Institute of Electronic Technology, Nanjing 210039, China 2.Department of Power Technology, Nanjing Lukou International Airport, Nanjing 211113, China MAXiande, XIEXiaofeng. Improved algorithm of discrete curvatures estimation for triangular meshes. Computer Engineering and Applications, 2014, 50（13）：196-200. Abstract：The Meyer method of calculating discrete curvatures for triangular meshes is concise in geometry and requires relatively small amount of calculation, but there is still potential to improve its results. The concepts of mean curvature structural vector and Gauss curvature structural angle are suggested, along with their geometric significance. Based on these analyses, the improved algorithms of the Meyer method are constructed, in which the main steps are simplified and the unnecessary errors are avoided. Simulation results show that the proposed algorithms effectively improve the calculation precision and efficiency of discrete curvatures estimation for triangular meshes. Key words：discrete curvature; triangulated grid model; mean curvature structural vector; Gauss curvature structural angle 摘 要：计算三角网格离散曲面曲率的 Meyer 方法几何意义简明，计算量较小，但其计算效果仍有进一步提高的潜 力。通过对 Meyer 方法的深入分析，提出了平均曲率构造向量和 Gauss 曲率构造角的概念，并指出了它们的几何意 义，在此基础上构造了对 Meyer 方法的改进算法。经分析，提出的改进算法精简了各个主要计算步骤，避免了不必 要的计算误差。仿真计算结果表明，改进算法是有效的，提高了三角网格离散曲率的计算精度和计算效率。 关键词：离散曲率；三角网格；平均曲率构造向量；Gauss 曲率构造角 文献标志码：A 中图分类号：TP391 doi：10.3778/j.issn.1002-8331.1207-0433 1 引言 网格的几何信息计算待求点处的曲率。 曲率是曲面上一点附近的重要局部几何性质，对于 有精确解析形式的曲面模型，其上任一点处的曲率可以 根据其坐标，由经典的微分几何方法的解析公式直接 计算[1]。对于没有精确解析形式，而只有表面离散点坐标 数据及离散点之间拓扑关系的曲面模型，则只能完全借助 于待求曲率的点附近各点的分布以及相互之间的拓扑 关系，采用数值方法求解。求解曲面上一点处的曲率是 获取该曲面几何信息的重要步骤，在离散几何、计算几 何、计算机图形学、虚拟现实技术等许多领域广泛应用。 由离散曲面上一点的邻近局部点计算离散曲率，一 般先把离散点划分为网格曲面模型，根据待求点各邻近 和 k Chen 等[2]在 1992 年提出，将经典微分几何和曲面论 之间 中描述法曲率 n、主曲率 k 1 关系的 Euler 公式改写为参数形式，并以圆弧拟合法和 最小二乘法求解法曲率参数，从而计算主曲率。此方法 容易实现，但圆弧拟合对真实的离散曲面而言过于粗 糙，精度不高。 与主方向 e 1 2 和 e 2 Taubin 等[3]在 1995 年根据曲面微分几何特性，以 n、 积分构造一个矩阵 B ，并证明 B 的特征向量就 和 e 。通过待求曲率 和 e e 1 2 是 n、e 1 的顶点的邻近三角面加权求和来估算矩阵 B ，进而计 算主曲率。 ，特征值中就包含了 k 、k 1 2 2 作者简介：马献德（1983—），男，博士，工程师，研究领域为数字信号处理，自动控制；谢小峰（1975—），男，工程师，主要研究领域 为自动控制。E-mail：mxdmxd@gmail.com 收稿日期：2012-07-29 修回日期：2012-11-26 文章编号：1002-8331（2014）13-0196-05 CNKI 网络优先出版：2013-01-15, http://www.cnki.net/kcms/detail/11.2127.TP.20130115.1140.006.html 

马献德，谢小峰：三角网格曲面上离散曲率改进算法 2014，50（13） 197 Meyer 等 [4]在 2003 年 采 用 离 散 化 Laplace-Beltrami ，采用离散化 Gauss- h ，并进而求解主曲率。 算子并积分来直接计算平均曲率 k Bonnet定理直接计算 Gauss 曲率 k g 处的一个法曲率可由点 p 1 Dong 等 [5]在 2005 年提出基于弧长参数微分计算法 及其附 的法向量及两点各自坐标估算。在此基础上 曲率的方法，即点 p 1 近一点 p 借助上述 Chen 等提出的方法求解主曲率和主方向，克 服了 Chen 法中圆弧拟合有时导致很大误差的问题。 2 此外，还有 Moreton 等 [6]提出的拟合法，Chen 等 [7]提 出的重心加权法，Page 等[8]提出的法向量投票法，柯映林 等[9]提出的通过 Shepard 曲面插值点线性组合法等。 其中，Meyer 的方法具有简明的几何意义，简单易 行，计算结果能达到一定的精度，且计算量较小。文献 [10]对多种算法进行了对比分析，认为 Meyer 方法的计 算效果最优。 经分析发现，此方法的主要计算步骤可以通过转换 得以简化，在计算精度和效率上可得到改进。 2 Meyer 方法概述 Meyer 提出的计算平均曲率 k 和高斯曲率 k 的方 g h 法[4]为： ( p k h ( p k g i i ) = 4A 1 Mixed ) = A 1 Mixed å Î N ( p 1 é ëê p 2π - å ) j i (cot α + cot β )v ij ij · n ù ij ûú （1） （2） θ æ èç 的所有直接邻近点（与 p Î N ( p 1 ö ø÷ p ) j j i i i 点处的曲面单位法向量。 v 共边的点） = p - p i j ij ， ) 是 p ( p 其中，N 1 的集合，n 为 p i i 其余依此类推。符号 · 表示向量内积运算。各参数的 定义如图 1 所示。 p i θ j β ij p j α ij 图 1 离散曲率计算参数示意图 A A Mixed Mixed 称为混合面积，其定义如下： p i å A( p = p j ) k （3） A( p p i p j k ) =  π/2||Ðp p j p k i  π/2 （4） p j i Î N ( p 1 ì ïï í ïï î )p ADijk ADijk ADijkV Î N ( p 1 k i k )p Î N ( p 1 /2Ðp p i j j )  π/2 p k p /4Ðp p i k j otherwise 是 Dp p p 的面积，而 A i ∆ijk 为锐角三角形时，图 2 中四边形 oh k j ∆ijkV 为 Dp p p 当 i j h p i k k j 的面积 其中，A Dp （点 o 为 Dp p p k i j p j p k i 的外心）。 h j O p k p i h k p j 图 2 锐角三角形的 Voronoi 面积 (A ) 示意图 ∆ijkV 上述公式需要多次计算三角函数值、反三角函数 值，由于涉及到三角形的边长、面积的计算，还需要多次 开平方，计算中容易引入截断误差和舍入误差，且计算 效率较低。以下试图精简计算步骤，提高计算精度和计 算效率。 i 假设 p 共有 m 个直接邻近点，从而有 m 个邻近三 的三角面），p 及其所有直接 角面（即有一个顶点为 p i 处 的 曲 面 单 位 法 向 量 n 均 已 邻 近 点 的 直 角 坐 标 ，p 知。研究目标是根据这些已知条件，以最精简的计算步 骤 完 成 以 Meyer 方 法 计 算 p 的 离 散 平 均 曲 率 和 离 散 Gauss 曲率的过程。设 + cot β )v i i i （5） ij ij （6） （7） j = (cot α v v = å = å Î N ( p 1 θ p j Σ ij v j ) i p Î N ( p 1 i ) j θ j 把 v 称为平均曲率构造向量，把 θ 角。可见 v、θ 称为 Gauss 曲率构造 的计算是算法中的主要步骤，决 Σ 、A Σ ∆ijkV 定了整个计算过程的计算量。以下分别寻找这三个参 数的改进算法。 3 对平均曲率构造向量的改进算法设计 3.1 平均曲率构造向量的普通算法简述 一般通过计算三角形三边长来计算公式（5）中的两 个余切值，则计算 v （1）计算 a = |v j | ，b = |v ij 的步骤如下： | ，c = |v kj ik | ，s = (a + b + c)/2 。 （2）采用三角形的性质来计算三角形内角的余切 值。对于 cot β ij 也按同样方法计算： s 1 cot α = s(s - a) ，s - s = (s 1 ij 2 = (s - b)(s - c) ，s 2 )s 3 = (s 1 s 2 )1/2 3 （3）采用公式（5）计算 v 。 j 以上步骤共 17m 次实数乘法，2m 次实数除法，4m 次实数开平方（加减法计算量很小，实数乘除 2 的整数 次幂可由位运算实现，故不考虑其计算量）。 3.2 平均曲率构造向量几何意义 首先分析平均曲率构造向量 v 的几何意义。把平 

198 2014，50（13） Computer Engineering and Applications 计算机工程与应用 均曲率计算公式中涉及的各个角度变换到一个三角网 格内，如图 3 所示。 p k α ij p i θ jk β ki p j 图 3 Meyer 法中各角度和向量变换示意图 比较图 1 和图 3 可知，v 可改写为： v = + v cot α (v å ij ij cot β ki ik ) （8） p Î N ( p 1 i j )p Î N ( p 1 i k )p Î N ( p 1 j ) k 设 v j = v ij cot α + v ik ij cot β ki 则有 v = å v ' j p i p k α ij β ki h p j h′ 图 4 v ' 的几何意义示意图 j 3.3 平均曲率构造向量的改进算法分析 由公式（9）可得： •v )2 ij |v ´ v | = |v |2|v |2 - (v ij ij ik ik ik 因而得改进算法步骤： （1）计算 u •v ，u = v = v •v ，u i ik a2 = |v ij |2 = v ij ij k •v ij jk ik ，b2 = |v = v j |2 = v •v kj •v ik ij ik ik 每 2 个相邻三角形有 1 条公共边，因而在 N 1 ( p i ) 内 a2 和 b2 的总计算量为 3m 次实数乘法。 （2）计算 l = |n ijk '| = 2A λ = 1/(2A ∆ijk ∆ijk ) = 1/l = |v ´ v | = (a2b2 - u2 i ij )1/2 ik 该步骤实际上附带得到了 Dp p j p k i 的面积，可在后 续的计算过程中使用。 （3）计算 w = λu k ，v ' = w i k v ij + (λu )v j ik k 以上步骤共有 22m 次实数乘法，m 次实数除法，m 次实数开平方。 4 对 Gauss 曲率构造角的改进算法设计 4.1 Gauss 曲率构造角的普通算法简述 通常，Gauss 曲率构造角按如下步骤计算： （1）对 p 的每个直接邻近三角形（即有一个顶点为 i 的三角形），根据下式计算角 θ ： p i θ = arccos v |v ik ik •v ||v ij ij | = arccos v •v ik ab ij 其中 a 和 b 在平均曲率构造向量的计算过程中已算出， 不必重复计算。 （2）求和： = å θ θ Σ j p Î N ( p 1 i ) j 以上步骤共有 4m 次实数乘法，m 次实数除法，m 次实数反余弦运算。 4.2 Gauss 曲率构造角的改进算法分析 利用了正切函数及两角和的正切公式，所得的改进 算法步骤如下： （1）对于 p i 的每个直接邻近三角形，计算 cot θ 的 j 值（或 θ j = π/2）。根据公式（10）有： j p )p Î N ( p 1 注意 v ' 与 v j k i j i )p Î N ( p 1 Î N ( p 1 并不相同，但对于同一个 v ，求和号中 ) k j v ' 的数目与 v j 的数目是相等的，均为 m 。因为 j v |v kj ki •v ||v | ，sin α = ij |v kj |v ´ v ||v ki | | （9） ki kj 表示向量 v 与 v ki kj ki kj 的外积运算，并依此类推。 cos α sin α ij = ij v |v •v kj ki ´ v kj ki | （10） cos α = ij 其中，v ´ v ki kj 所以 cot α = ij 同理 cot β = ki v |v •v ji jk ´ v ji | 注意到 |v ´ v kj •v jk (v ik v ' = j jk | = |v ki + (v ´ v ij | ij ji •v ij )v |v ik ´ v p j i p k ´ v )v kj jk ik | = |v ´ v ik ij | = 2A ∆ijk ，可得 = v kj ik ´(v 2ADp ) ij = v kj ) ij ´(v |v ik ´ v | ik ´ v 如图 4 所示，在平面 p p j p k i ij 内过点 p （11） 作直线 p p k j i 的垂线 p i h ，垂足为 h 。延长线段 p h 到 h′ ，使得 i |p i - h′| = |v | kj 则有 v ' = h′ - p j 由此可见，v 实际上是 p i 点的所有邻近三角面片 i 点所在三角形的对边的正交向量的加权和，权值 点对边的长度。这就是平均曲率构造向量 v 的 中，p i 就是 p 几何意义。 i 

马献德，谢小峰：三角网格曲面上离散曲率改进算法 2014，50（13） 199 cot θ = j v |v •v ik ij ´ v ik ij = λ(v •v ij ik ) = λu i | （12） （2）计算 cotθ Σ 的值（或得到 θ Σ = (2r + 1)π/2 ，r 为非 负整数），依次计算： cot æ ç è t + 1 ö ÷å θ ø k = 1 k = cot æ ç è t å k = 1 + θ θ k t + 1 = ö ÷ ø t cotå θ k = 1 t cotå k = 1 cot θ k t + 1 - 1 + cot θ θ k t + 1 其中 t 的取值依次为 12m 。注意须先处理分母为 零的特殊情况，同时根据参与计算的两个余切值和所得 结果的正负号，确定所得结果对应的角度所在的区间。 最后得到 cot θ = x 。 Σ （3）计算 θ Σ = arccot(x) + m π 0 0 其中，m 确定。一般而言，m 次加法实现。 0 为非负整数，由上述得到的 θ 所在的区间而 Σ 的值不大，此处的 m π 可用少数几 0 以上步骤共有 2m 次实数乘法，m 次实数除法，1 次 实数反余切运算。 5 对混合面积的改进算法设计 5.1 混合面积普通算法简述 A ∆ijk 一般可采用如下方法计算： ADijk A ∆ijkV = s(s - a)(s - b)(s - c) = s 1 s 2 的计算步骤如下： （1）设 Dp p p 外接圆半径为 R ，则： i k /16 j R2 = a2b2c2s2 3 （2）计算图 2 中线段 oh | oh = R2 - b2/4 ， | | oh k ： （3）计算面积 A | j j A ∆ijkV = (a|oh j ∆ijkV | + b|oh |)/4 k 和 oh k 的长度。 = R2 - a2/4 注意到在平均曲率构造向量的普通算法中，实际上 已经计算出 a2、b2、c2 和 s2 。因而以上步骤共有 5m 次 3 实数乘法，2m 次实数开平方（假设所有三角面都是锐角 三角形，下同）。 5.2 混合面积改进算法分析 如 图 5 所 示 ，连 结 p o 并 延 长 到 点 p ，使 得 m p p k j ||oh ，连结 p p 与 p 。推导可得 m (v ik ADmjk = j ) •v ij kj k m •v )(v 4ADp jk p p k j i m ì ï ï í ï ï î i p = 1 2 ) = 1 8 ADijkV = 1 4 (ADijk + ADmjk (l n ' ijk + w u ) j k 故共需 m 次实数乘法。 λu u j k （13） p i θ h k β p j h j p k α O p m p n 图 5 Voronoi 面积改进算法构造示意图 5.3 混合面积计算方法的判断 前文已提到，计算混合面积时，对不同的三角形，其 计算方法不同。根据前述计算，有 A( p p i p j k ) = ì ïï í ïï î ADijk ADijk ADijkV /2u /4u i  0  0||u j k otherwise  0 （14） 其中 u 、u i j 和 u k 已在前述计算步骤中得到，因此对每个 三角面片，决定采用哪种方法计算其混合面积时，并不 需要额外的计算量。 综上所述，当所求点处的直接临近三角形全为锐角 三角形时，普通算法共需 26m 次实数乘法，3m 次实数 除法，6m 次实数开平方，m 次实数反三角函数；而改进 算法共需 25m 次实数乘法，2m 次实数除法，m 次实数 开平方，1 次实数反三角函数。 可见，改进算法中除了乘法次数与普通算法基本持 平外，计算效率很低、耗时很长的除法、开平方和反三角 函数的次数减少，因此提高了计算效率，同时更重要的 是由于减少了计算量，避免了不必要的截断误差和舍入 误差，从而有望减小计算误差，提高计算精度。 6 改进算法的效果 为验证所提算法的有效性，构造离散的椭球面、椭 圆柱面、椭圆抛物面、双曲抛物面、双叶双曲面、单叶双 曲面等 6 种离散标准二次曲面，分别按普通算法和所提 改进算法计算对比。步骤如下： （1）给定曲面形式、几何特征和空间范围，随机生成 满足这些条件的离散点坐标及其三角网格拓扑结构，并 计算各点处平均曲率和 Gauss 曲率的真实值。 （2）计算各点处的法向量真实值，在上述算法中的 离散法向量均直接采用真实值，这样可避免采用离散算 法计算法向量带来的误差。 （3）在同等计算条件下，分别采用普通方法和所提 改进方法来实现 Meyer 算法，得到两种方法的相对误差 （剔除因真实曲率为 0 而无法计算相对误差的点）；并记 录两种算法所耗费的计算时间，得到计算精度和计算效 率的评价。 采 用 的 计 算 软 硬 件 平 台 为 ：Intel CoreTM 2 Duo T5870@2.00 GHz CPU，4.00 GB 内存，Windows 7 操作 ® 

200 2014，50（13） Computer Engineering and Applications 计算机工程与应用 系统，Visual C++ 2008 开发平台。 7 结束语 曲率相对误差的统计直方图如图 6 所示。普通方 法的平均曲率相对误差的均方根值为 7.722%，Gauss 曲 率相对误差的均方根值为 7.712%；改进方法平均曲率 相对误差的均方根值为 1.066%，Gauss 曲率相对误差的 均方根值为 1.058%。 2.5 2.0 1.5 1.0 0.5 4 0 1 / 量 数 占 所 2.5 2.0 1.5 1.0 0.5 5 0 1 / 量 数 占 所 0 －0.5 0 曲率相对误差 0.5 0 －0.2 －0.1 0 0.1 曲率相对误差 0.2 （a）平均曲率（普通方法） （b）平均曲率（改进方法） 3 2 1 4 0 1 / 量 数 占 所 2.0 1.5 1.0 0.5 5 0 1 / 量 数 占 所 0 －0.5 0 0.5 曲率相对误差 1.0 0 －0.2 －0.1 0 0.1 曲率相对误差 0.2 （c）Gauss 曲率（普通方法） （d）Gauss 曲率（改进方法） 图 6 计算精度对比（相对误差直方图） 对一个离散点的平均曲率与 Gauss 曲率的计算耗 时的统计直方图，如图 7 所示。普通方法的耗时均方根 值为 1.222 ms，改进方法耗时均方根值为 1.036 ms。 4 0 1 / 量 数 占 所 15 10 5 0 4 0 1 / 量 数 占 所 15 10 5 0 5 10 计算耗时/ms （a）普通方法 5 10 计算耗时/ms （b）改进方法 图 7 计算效率对比（计算耗时直方图） 由 计 算 结 果 可 见 ，所 提 的 改 进 方 法 大 大 提 高 了 Meyer 算法对离散平均曲率和 Gauss 曲率的计算精度， 同时还提高了其计算效率。 通过分析 Meyer 所提的离散曲面一点处平均曲率和 Gauss 曲率的算法，提出了平均曲率构造向量和 Gauss 曲率构造角的概念，分析了其几何意义，并据此提出了 提高计算精度和计算效率的改进算法。对 6 种典型离 散标准二次曲面的仿真计算结果，表明了本文改进算法 可大大提高曲率计算精度，同时还可提高计算效率。 参考文献： [1] do Carmo M.Differential geometry of curves and surfaces[M]. Englewood Cliffs，NJ：Prentice-Hall，1976. [2] Chen Xin，Schmitt，F.Intrinsic surface properties from sur- face triangulation[C]//Proceedings of the European Con- ference on Computer Version.Heidelberg：Springer-Verlag， 1992，588：739-743. [3] Taubin G.Estimating the tensor of curvature of a surface from a polyhedral approximation[C]//Proceedings of the Fifth International Conference on Computer Vision.Pisca- taway，NJ：IEEE，1995：902-907. ö [4] Meyer M，Desbrun M，Schr der P，et al.Discrete differential triangulated 2-Manifolds[J].Visu- geometry operators for alization and Mathematics，2003，3：35-57. [5] Dong Chenshi，Wang Guozhao.Curvatures estimation on triangular mesh[J].Journal of Zhejiang University Science， 2005，6A：128-136. [6] Moreton H P，Sequin C H.Functional optimization for fair surface design[C]//Computer Graphics Proceedings，Annual Conference Series（ACM SIGGRAPH），Chicago，Illinois， 1992：167-176. [7] Chen S，Wu J.Estimating normal vectors and curvatures by centroid weights[J].Computer Aided Geometric Design， 2004，21（5）：447-458. [8] Page D L，Sun Y，Koschan A F，et al.Normal vector voting：crease detection and curvature estimation on large， noisy meshes[J].Graphical Models，2002，64（3/4）：199-229. [9] 柯映林，陈曦.基于 4D Shepard 曲面的点云曲率估算[J].浙 江大学学报：工学版，2005，39（6）：761-764. [10] 方惠兰，王国瑾.三角网格曲面上离散曲率估算方法的比 较与分析[J].计算机辅助设计与图形学学报，2005，17（11）： 2500-2507.