基于离散指数映射的网格曲面上曲线设计与复用

韩 林 刘 斌

华侨大学，厦门，361021

0 引言

在几何造型中，曲面上的曲线设计是一项重要的内容。参数曲面上的曲线表达研究越来越多［1-3］，但对于网格曲面上的曲线，尚没有理想的表示和设计方法。随着计算机辅助设计和现代制造业的发展，由平面和参数曲面表示的传统规则特征已远远不能满足现代工业产品设计在外形上的要求，尤其是在计算机动画、数字影视制作，以及柔性制品的交互设计中，参数曲面造型更是存在较多空白。三角网格模型能够表示任意复杂形状，以其独特优势成为主流表示方法之一，故研究网格曲面上曲线的设计方法与配套技术具有适用价值和理论意义。

国内外学者对曲面上曲线的生成方法已有一些研究：Hofer等［4］给出了一种光滑流形曲面上三次样条曲线的生成方法；Rodriguez等［5］、Park等［6］将De Casteljau算法应用于简单的流形曲面；Morera等［7］给出了测地Bezier曲线的生成方法；杜宏云［8］在文献［6］的基础上给出了测地自由曲线的生成算法；朱文明等［9］从细分的角度构造出流形曲面上的细分曲线。刘斌等［10］提出了一种插值于流形网格曲面上给定点列的测地B样条曲线生成方法，该方法以网格曲面上两点间的最短测地线代替欧氏空间中的两点间直线，将欧氏空间中的德布尔算法拓展到曲面空间，形成与经典B样条曲线具有统一表示形式的曲面上自由曲线的表达形式。由于测地B样条具有局部控制性、凸包性等优点，故本文以文献［10］所提出的方法来生成网格上的自由曲线，进行交互操作。

欲实现测地B样条曲线在弯曲空间中的基本交互操作，使其类似于B样条曲线在欧氏空间中的基本操作，例如平移、旋转、缩放等，有效方法之一就是网格参数化的方法。三角网格参数化是对这些三角网格几何和拓扑信息进一步处理的基础［11］，各种各样的参数化方法研究一直是计算机图形学的热点。通过分析各种参数化方法的优缺点，发现离散指数映射近似法［12］计算高效，而且在局部范围内保留了三维空间中每一点的距离与方向不变的特性，将其应用于实时交互操作与复用设计，具有独特的优势。

本文在分析上述领域研究现状的基础上，把指数映射参数化的方法引入到网格曲面上自由曲线的设计中，提出一种网格曲面上测地B样条曲线的复用操作方法（主要包括平移、旋转、缩放等），为优秀设计成果的保存和重用提供了支持，提高了形状设计的智能程度和效率。

1 离散指数映射

1.1 指数映射

在曲面S：r＝r（u1，u2）上取定一点p，p点的切平面为Tp，任取切向量v∈Tp，作测地射线Cp，v从p点出发并且以v／｜v｜为初始切向，则Cp，v由p和v／｜v｜唯一确定，如图1所示。

图1 曲面S上点P处的指数映射

取Cp，v的正向弧长s参数化ui（s），i＝1，2，使p点在S 上的曲线坐标为（u1（0），u2（1））。定义映射：

即像点Qs＝r（u1（s），u2（s））是CP，v上从p点出发且经过弧长s后到达的点。由此定义映射为

则此映射称为曲面S上点p处的指数映射［13］。

指数映射expp本质上就是在p点附近的曲面上产生一个二维的坐标系统，将曲面S上的点映射到p点的切平面Tp上，并称p点为种子点。对任一单位向量v∈Tp来说，都存在一个用弧长来参数化的测地线gv，使得gv（0）＝p 且g′（0）＝v。对p点周围相邻的任一邻域点Q而言，都会有唯一的测地线gv通过它。因此，Q可用极坐标（ρ，φ）映射到Tp上，ρ是从p到Q的测地线距离，φ是v在TP上的极角，这就是测地极坐标表示形式。极坐标可在Tp上被任意一个正交基底｛η1，η2｝表达成法坐标（u，v），也可以将其看成二维向量e＝（u，v），此二维向量称为测地线向量。本文中，这2个概念将不加区分地使用，法坐标为坐标表示形式，测地线向量为向量表示形式。指数映射参数化的关键步骤就是求取曲面上任意一点到种子点的测地线向量。

1.2 离散型指数映射逼近

对于网格曲面，为了近似得到曲面S上点p处的指数映射，若按定义，需要计算曲面上其他点q相对于点p的测地线距离和极角。测地线的计算成本较高，因此，本文采用Schmidt等［12］提出的离散指数映射逼近法进行测地线的计算，该方法并不具体计算点q到点p的曲面测地线和极角，而是先由Dijkstra算法产生分段线性测地线向量，再通过这些测地线向量叠加，计算点q在点p切平面上的法坐标。

假设网格曲面上有3点，分别为p、r、q，如图2所示。p到r的测地线向量up，r和r到q的测地线向量ur，q已知，而p到q的测地线向量未知。目标是求取up，q＝expp（q），即点q在切平面Tp上的法坐标。在线性系统中存在

图2 离散指数映射逼近

由于（up，q-up，r）未知，用ur，q的一个相关量来代替（up，q-up，r）。具体方法是，将原来用Tr的三维基底（xr，yr，zr）表示的ur，q用Tp的三维基底（xp，yp，zp）来表示。首先，设点p、r处的法矢分别为np与nr，将Tr沿着nr与np的叉积方向（nr×np）旋转两矢量之间的夹角α，此时Tr与Tp共面，Tr的基底变成（x′r，y′r，np）；然后，再将Tr沿着np的方向旋转角度θp，r＝arccos（x′r·xp），从而Tr与Tp具有相同的三维基底。经过两次旋转后，Tr上的ur，q就可以用Tp的基底来表示。由于up，r与ur，q都是二维向量，可以忽略第一次旋转对ur，q造成的影响，因此，up，q的近似值可表示为

式中，Rot2 D（θp，r）表示二维旋转变换。

虽然用式（2）通过测地线向量旋转、叠加求取的up，q与实际的三维网格上测地线向量的值会存在误差（图2），但是若不谋求网格全局参数化，只考虑局部映射，精度是能够保证的，原因是此方法采用了式（2）这样的矢量叠加，而不是以标量运算。通过式（2）的延伸，就可以近似求出点p局部邻域网格上每一个顶点在Tp上的测地线向量。

1.3 法坐标归一化

基于指数映射的曲线交互设计方法需要确定种子点p和区域半径R，如图3所示。种子点作为算法向外扩张的起始点，种子点的切平面与参数平面共面；区域半径R是种子点到区域范围内最远点的测地线距离，R既决定了区域的范围，又是算法的停止条件。计算出区域内的所有顶点的测地线向量后，便得到了区域内顶点的法坐标，然后将法坐标进行几何缩放和平移变换。其中，比例变换矩阵为

平移变换矩阵为

通过2次变换将所有顶点的法坐标映射到［0，1］×［0，1］的参数空间中。

图3 法坐标归一化

2 测地B样条曲线交互操作

用户在网格曲面上进行测地B样条曲线的编辑时，希望看到测地B样条曲线在网格曲面上的滑动、旋转和缩放等变换，且要满足实时交互性。基于这些要求，本文首先在网格模型上输入型值点，然后插值测地B样条曲线（图4a）；其次，在网格模型上选出一块区域，该区域（称之为源区域）必须完全覆盖测地B样条曲线，按照指数映射参数化的方法计算该区域内测地B样条曲线控制顶点的法坐标（图4b）；最后，通过选择或改变网格模型上区域（称之为兴趣区域）的位置、方向或大小实现测地B样条曲线的重用编辑（图4c）。结果如图4d所示。

图4 曲线交互操作复用流程

2.1 测地B样条曲线的生成

在网格曲面上生成k次测地B样条曲线的算法思路是：首先给定三角网格曲面模型G和位于网格曲面上的m＋1个型值点Qj∈G（j＝0，1，…，m）；其次，由各型值点之间的最短测地距离按照积累弦长参数化方法进行数据参数化，得到各型值点的参数值，类似于经典B样条插值理论，在端点处取k＋1重节点的固支条件，且取规范定义域，这样可得欲求曲线的节点矢量，然后由插值条件确定欧氏空间中各控制顶点；最后，将得到的控制顶点投影到网格曲面上，按扩展德布尔算法得到初始测地B样条曲线，计算各插值点到其在曲线上对应点的测地距离，即插值误差，按照误差反向补偿策略，更新型值点，重复上述过程，直至误差满足精度要求为止。生成B样条曲线的同时，把计算出的控制顶点按顺序存取，作为后续步骤的重要信息。

2.2 源区域选择和控制顶点参数化

2.2.1 源区域参数化

生成测地B样条曲线之后，需要进行源区域的选择及参数化源区域，其关键步骤是求取源区域内任意一点到种子点的测地线向量。

3D网格上的每一个顶点与其所有1-ring邻域顶点之间的测地线距离，即为连接两点的边的长度，因此其所有1-ring邻域顶点的测地线向量可以很轻易地计算出来。如图5所示，对于顶点p的1-ring邻域顶点q，利用旋转的方式将向量vp，q沿着vp，q与p的切平面法向量np的外积方向旋转至p的切平面上，这样既可以得到q点的方向，又能保持p、q之间的距离不变。然后将旋转过后的三维向量，以切平面上一组正交基底｛xp，yp｝表示成二维向量up，q，up，q即为q点在p的切平面上的测地线向量。

图5 网格顶点1－ring邻域的测地线向量

本文以每个顶点的任意一个1-ring邻域顶点来建立基底。假设一平面Tp由一组正交基底｛η1，η2｝所生成，则在Tp上的某点q与平面的原点p所形成的向量l可表示为

式（3）中，因u、v为标量，η1、η2都是单位向量，故u、v可分别表示为

二维向量up，q＝（u，v）即为q点在Tp上的测地线向量。用同样的方法可求取源区域内所有顶点1-ring邻域的测地线向量。

按照前述离散指数映射理论，若种子点到源区域任意顶点的路径已知，便可利用向量的旋转叠加得到种子点到源区域任意顶点的测地线向量。

如图6所示，假设种子点为p，利用Dijkstra最短路径方法得到p到网格上源区域内任一顶点q的最短路径｛p0，p1，…，pm｝，p0＝｜p｜且pm＝｜q｜，之后沿着最短路径不断地将每个切平面上的测地线向量转换到Tp的基底上，且其q点的测地线向量可表示为

图6 种子点p到网格上任意顶点q的测地线向量

须注意的是：二维向量upi，pi＋1是用直接转换成起始点的切平面Tp的基底来表示的，而非前一个切平面Tpi-1的基底。如果通过前一个切平面进行转换的话，则会增大误差积累［12］。

计算出源区域内种子点到所有顶点的测地线向量后，即可对其进行法坐标归一化处理，使所有顶点的法坐标映射到u∈ （0，1），v∈ （0，1）的范围中。网格曲面上三角面内部点的法坐标可通过其所在的三角面片的重心坐标求取。

2.2.2 曲线控制顶点参数化

曲线控制顶点一般不会刚好位于网格曲面的顶点上，一般情况下都位于三角面片内部（包括位于三角面边上），这时可用其所在三角面片3个顶点的法坐标和该点的三角形重心坐标求取该点的法坐标。如图7所示，设曲线控制顶点P位于网格曲面的三角面片△ABC内，3个顶点的测地向量分别为a、b、c，P点在三角面片内，A、B、C3个顶点的重心坐标u、v、w可分别表示为

则P点的测地线向量可以表示为

计算曲线所有控制顶点的法坐标（测地向量）并将其存储备为曲线操作的基础数据。

2.3 参数匹配

无论是曲线的滑动迁移，还是旋转缩放，在操作过程中，测地B样条曲线各控制顶点的法坐标应始终保持不变，改变的只是兴趣区域的位置、方向和大小。选定兴趣区域后的主要任务就是从兴趣区域中找到源区域的对应点，使其与源区域一一对应。

图7 △ABC的内点P的重心坐标

图8 源区域G1与兴趣区域G2的参数化

如图8所示，假设源区域为G1，兴趣区域为G2，B样条曲线上n个控制顶点为pi∈ ｛p0，p1，…，pn｝，两区域的方向和大小都没有改变，只是位置平移。两区域是同一网格的不同位置，很难保证两区域的网格具有相同的三角面片个数和相同的连接关系，为此提出一种拓扑不同的三角网格区域参数匹配的快速算法，具体步骤如下：

（1）按照参数化源区域G1的方法，将兴趣区域G2参数化。保存G2内所有顶点的法坐标，以此建立k-d树。搜索G2参数域内与pi的法坐标（ui，vi）最近的顶点，以该顶点的1-ring邻接三角形为候选搜索集，求出pi在兴趣区域G2的对应点p′i所在的三角面片 △abc。

（2）取出△abc 3个顶点的法坐标，用这3个法坐标建立虚拟三角形，称之为该三角面片的参数三角形。由第（1）步可知，（ui，vi）一定落在△abc的参数三角形内部（包括落在三角形边上的情况），利用类似于前述重心坐标的求解方法，得到（ui，vi）在参数三角形中的重心坐标。

（3）由重心坐标和三角形3个顶点的三维空间坐标，即可计算出p′i的三维空间位置，即pi在兴趣区域中的对应点。

计算出B样条曲线各个控制顶点在兴趣区域中的对应点之后，按照拓展德布尔算法再生成测地B样条曲线，就完成了曲线的滑动平移操作。

每一个兴趣区域都是基于种子点、方向及半径产生的，只要改变兴趣区域的半径范围，再以上述步骤重新计算B样条曲线的控制顶点的法坐标，就可以实现缩放的效果。测地B样条曲线的旋转则是将兴趣区域种子点切平面的正交基底｛η1，η2｝沿着切平面的法向量旋转一个角度，再以上述步骤重新计算B样条曲线的控制顶点的法坐标的过程。

3 试验结果与分析

本文通过实例展示算法的结果，主要考虑了算法在各种曲面上的适应性和视觉效果，并对曲线的复用交互设计速度进行了探讨分析。

3.1 算法的适应能力及其视觉效果

本文算法使用Microsoft Visual Studio 2008开发工具和OpenGL函数库，在CPU为2.26GHz、内存为2.0GB的PC机上实现，部分结果如图9所示。图9a中的曲线位于半球面上，用于测试曲线变换后的视觉效果；图9b中的帽子模型中部凸起，用于测试曲线在曲率变化较大的曲面上的适应性；图9c和图9d为圆鼓和圆台模型，用于测试各种曲线在具有尖锐特征的曲面上的适应性。试验结果表明，本文提出的曲面上测地B样条曲线的复用方法视觉效果良好，曲线对曲面的适应性较强。

3.2 时间效率

良好的视觉效果和较强的适应性是算法的基本要求，要应用于实际工程，还须满足实时交互的要求，即算法时间应该控制在一定范围内。为此，选择图10所示的Shoe模型和Venus模型作为试验对象，试验不同曲线图形在两模型中的重用时间。另外，针对本文提出的仅计算局部区域的快速参数匹配算法，与计算整个网格区域的非快速算法进行了时间效率比较，试验数据如表1所示。

从表1中的数据可以看出，快速参数匹配的时间要比非快速参数匹配的时间低一个数量级，网格规模愈大，快速参数匹配的优越性愈明显。图10a图形和图10b图形都位于Shoe模型上，但图10b图形的匹配时间为图10a图形的2倍，因图10b的控制顶点为图10a的2倍，变换时须增加匹配次数，故时间变长。匹配时间与控制顶点的个数之间基本上呈线性关系，通过Venus模型上的图10c图形和图10d图形可以进一步给予证实。图10a图形和图10c图形都为10个控制顶点，但后者所在的参数区域的顶点个数是前者的5倍多，导致图10c图形的匹配时间线性增加。图10d图形的匹配时间比图10b长，也是此原因。

图9 各种网格曲面上测地B样条曲线图形的交互操作

4 结论

（1）将指数映射理论引入网格模型的局部参数化，实现网格曲面上曲线的迁移、旋转和缩放等操作，为曲面上曲线的复用和再设计奠定了基础。

图10 网格模型上测地B样条曲线图形的复用操作

表1 曲线插值迭代时间与模型的关系

（2）所采用的局部参数化方法计算效率高，变形小，不依赖于父网格模型的整体规模，能够满足实时交互设计要求。

（3）采用测地B样条作为网格曲面上的曲线建模工具，曲线对曲面的适应性强。

（4）离散指数映射计算速度快，但由于存在积累误差，不适宜于大范围网格区域的计算，如何减小积累误差的影响，有待于进一步的深入研究。

［1］蔺宏伟，王国瑾.光滑曲面上的G1插值曲线［J］.计算机辅助设计与图形学学报，2003，15（5）：541-546.

［2］冯结青，彭群生.基于插值的Bernstein多项式复合及其曲线曲面应用［J］.软件学报，2002，13（10）：2014-2020.

［3］卫炜，周来水，王志国.NURBS曲面上的曲线精确表达［J］.南京航空航天大学学报，2008，40（1）：85-88.

［4］Hofer M，Pottmann H.Energy-minimizing Splines in Manifolds［J］.ACM Trans.Graph.，2004，23（3）：284-293.

［5］Rodriguez R C，Leite F S，Jacubiak J.A New Geometric Algorithm to Generate Smooth Interpolating Curves on Riemannian Manifolds［J］.LMS Journal of Computation and Mathematics，2005，8：251-266.

［6］Park F C，Ravani B.Bezier Curves on Riemannian Manifolds and Lie Groups with Kinematic Applications［J］.ASME Journal of Mechanical Design，1995，117：36-40.

［7］Morera D M，Carvalho P C，Velho L.Modeling on Triangulations with Geodesic Curves［J］.The Visual Computer，2008，24（12）：1025-1037.

［8］杜宏云.测地自由曲线及其性质研究［D］.南京：航空航天大学，2008.

［9］朱文明，邓建松，陈发来.应用保角映射构造流形上的细分曲线［J］.计算机辅助设计与图形学学报，2007，19（1）：48-53.

［10］刘斌，黄常标，林俊义，等.流形网格曲面上测地B样条插值［J］.机械工程学报，2011，47（19）：136-142.

［11］彭群生，胡国飞.三角网格的参数化［J］.计算机辅助设计与图形学学报，2004，16（6）：731-739.

［12］Schmidt R，Grimm C，Wyvill B.Interactive Decal Compositing with Discrete Exponential Maps［J］.ACM Trans.Graph.，2006，25（3）：605-613.

［13］王幼宁.微分几何讲义［M］.北京：北京师范大学出版社，2007.