杨 峰 廖文和 戴 宁 董光雷 韩厚年 郭保苏 孙玉春

1.南京航空航天大学，南京，210016 2.北京大学口腔医学院，北京，100081

0 引言

近年来，随着制造技术的飞速发展与物质产品的日益丰富，如何将CAD/CAM等现代制造技术用于生命本身质量的提高，逐渐得到人们更多的关注。在生物医学领域，需要采用人工方法修复病损组织器官的案例日益增多，常见的需要修复的生物体有人工关节、颅颌面修复体、口腔修复体、整形外科植入体等。由于人体组织器官结构的复杂性和临床医学个性化需求的多样性，要保证修复体与病人病损部位的精确匹配，就必须进一步提高CAD/CAM的设计与制造精度。为适应上述发展和要求，本文提出一种基于累积曲率的离散刀轨环点云数据预处理优化算法。

1 生物修复体的数据采集、建模及制造

1.1 离散网格曲面与网格数及点云数的关系

修复体的制作主要包括三维数据采集、模型构建和加工制造三个步骤。数据采集中得到的是许多离散的点云数据，为了提高匹配精度，在模型构建阶段往往采用离散网格去逼近实体模型，而在加工制造阶段，为了减小因曲面造型产生的逼近误差，往往直接对网格曲面生成刀位轨迹。离散网格曲面是许多空间多边形面片对实体模型的线性逼近，多边形面片的疏密决定了逼近精度的高低，同时会影响精确刀位轨迹的质量和生成效率。离散网格曲面通常采用截面线法直接生成刀位轨迹［1－2］，该方法以曲面每一层的二维截面轮廓点云数据为基本元素，通过对其进行等距、连接、过渡等几何处理得到刀位轨迹。在模型构建时，为了更精准地表示实物模型的细节信息，总是希望多边形面片密度越大越好，这样在二维截面轮廓生成过程中，就会获取到大量的高精度点云数据。图1a中的牙齿模型由6312个三角面片构成，可以很好地反映实体模型的细节信息；图1b中的牙齿模型由631个三角面片组成，基本可以反映实体模型的整体信息，但部分细节信息已经失真；图1c为该模型某二维截面轮廓；图1d为图1c的细节放大图，其中点云数据分布非常密，很多点都是冗余的。然而，庞大的轮廓点云数据量会给轮廓曲线的等距、过渡等后续处理及刀位数据输出带来严峻的挑战，最终造成许多较短的、不规则的加工路径出现，导致加工过程中的频繁抬刀及刀具在轨迹突变处的振动，势必降低加工精度和效率，这一点在高速加工中尤其显著。图2所示为点云数据对等距的影响。图2a为原始轮廓数据的等距效果图，图2b为图2a的局部放大图，不难看出，过密的点云等距后会出现很多难处理的自交等不规则现象。可见，高密度的点云数据并不总是必要的，相反，稠密的点云数据不仅会影响后续相关几何处理的精度和效率，而且会造成加工中的许多问题。因此，计算机辅助制造模块在生成二维截面轮廓后，为提高精度，需要适当地从截面轮廓点云数据中抽取对模型加工有用的信息，并按照一定的原则对轮廓曲线进行光顺处理，才能生成较好的刀位轨迹，从而满足加工要求。

图1 网格数量对离散网格曲面的影响

图2 点云数据对等距的影响

1.2 离散网格曲面点云数简化方法

有关点云数据简化的方法很多，主要有：包围盒法、平均距离法、基于曲率的简化方法、三维栅格细分法。大部分简化方法都用于曲面造型中，而对于二维截面轮廓点云数据的简化，最常用的是基于曲率的简化方法。Teh等［3］用计算点的离散曲率来删除冗余点，直接以点到对边的距离h除以对边的长度L作为近似曲率。该方法可有效去除冗余点，但当原始曲线变化比较缓慢时，有可能会误删除一些有效的特征点而产生过大的累积误差或者错误。Hamann等［4］在简化平面图形和数字化模型的基础上，提出将图形中的曲线因子进行分段，通过计算局部区域的曲线曲率值来选取保留的数据。该方法的缺点是点云数据简化的效率受选取的点云数目制约，也和最终的误差阈值有关。Martin等［5］采用先把所得的数据点进行均匀网格划分，再从每个网格中提取样本点组成优化点云链的方法实现数据精简。该方法精简效率高，但对物体外形尺寸和曲率变化不敏感。平均距离精简法［6］容易丢失细节信息，一般用于扫描线形式的点云数据。

1.3 离散网格曲面曲线光顺方法

传统的曲线光顺方法大致分为整体光顺法和局部光顺法。Liu等［7］给出了不带约束的全局优化光顺法，而文献［8］给出了带约束的全局优化光顺法。全局优化光顺方法的优点是可以一次完成对曲线的光顺，缺点是一次光顺的点太多，计算量大，难以保证个别坏点的特殊修正。同简化方法类似，光顺处理大多基于CAD研究，因计算繁冗，故二维截面轮廓的光顺多采用局部光顺法。Liu等［9］提出了一种离散曲率（D－curvature）光顺算法，并开发了相应的自适应算法。该方法利用寻优函数的约束确定坏点的修正量，但因其没有限制修正量的阈值，故其曲率难以满足加工精度的要求。Farin等［10］用绘制曲线曲率图的办法，人为地判断曲线的光顺程度，并指出需要修改的点，该方法直观易行，但人工交互繁杂耗时，且受屏幕分辨率的影响较大。Renz［11］用直接计算点列的一阶和二阶差分方法人为地判断坏点，并对其进行光顺，该方法的不足是没有考虑点的均布问题。

2 离散曲率的相关技术定义

2.1 点云数据的离散曲率

平面曲线的曲率计算公式为

本文把曲率计算转化为三角形内边的计算，并提出通过直线方程判断曲率正负号的离散曲率计算方法。如圆逼近法求曲率，是通过相邻3个点Pi－1，Pi，Pi＋1的圆来近似Pi点的密切圆，如图3所示。圆的曲率Ki为

图3 点的离散曲率计算

因由离散点列组成的曲线有凹凸之分，所以点的离散曲率也有正负之分。由于式（2）中的分子部分是把点Pi的坐标直接代入对边Pi－1Pi＋1的方程所得，所以是有正负的，且其正负直接决定点Pi的离散曲率正负。

2.2 二阶中心差分、差商光顺

设给定型值点列｛Pi｝（i＝0，1，…，n），则Pi的二阶中心差分、差商分别为

式中，Li为由P1到Pi的累积弦长。

适用于二阶差分光顺和差商光顺的点云分布条件如下：

（1）若点云数据分布比较均匀，则采用二阶差分法光顺点云数据。对于型值点Pi，调整其位置，使得Δ2Pi＝0，则可得Pi点的新位置：

（2）若点云数据分布不均匀，则需要用差商法代替差分法才能达到较好的光顺效果。对于型值点Pi，调整其位置，使得f（Pi－1，Pi，Pi＋1）＝ 0，则可得Pi点的新位置：

2.3 基于光顺算子和简化算子的误差带控制

由于本文的研究对象是二维截面轮廓点云数据，而经过等距、连接和过渡等几何处理后的轮廓数据将直接作为刀位路径用于实际加工中，所以，二维轮廓点云数据优化的最大误差应控制在允许的加工误差范围内。为满足加工精度要求，提高加工效率，本文取Pi优化后的位置为原来位置和理想位置的组合，并提出图4所示的加工误差带容差δ概念。

优化后Pi点的新位置为

图4 加工误差带示意图

式中，α为简化算子，α∈ ［0，1］；β为光顺算子，β∈ ［0，1］。

要确保加工精度，就必须使

式（6）中的容差δ由实际加工中的精度要求确定。

在优化过程中，所有点必须在误差带控制范围内，即在图4的优化过程中，二维轮廓点云数据的误差应该控制在内外两条实线组成的封闭轮廓范围内。如果新生成的点P″i不能满足式（6），则修改光顺算子α，采用递归的方式，直到P″i满足误差带控制要求为止。

图5为二维截面轮廓点云数据优化的整体示意图，可以看出，优化过程主要分为两步，即先简化，后光顺。简化分为粗简化和精简化两步：如果粗简化满足条件，则直接进入光顺处理，如果不能满足，需要精简化后进入光顺处理。整个处理过程都在误差带控制范围内。

图5 二维截面轮廓点云数据优化示意图

3 基于累积曲率的二维截面轮廓点云数据的简化算法

理想的简化算法总是希望用最少的点云数据表达最必要的信息，但在实际处理过程中通常根据某种简化原则来进行。对于一个给定的点云数据，在不同的应用领域有多种简化准则，主要包括：保留数据点的个数比重、数据点间的空间距离、曲率、密度及重建误差等准则。对于二维截面轮廓点云数据，由于优化后的数据直接用于生成刀轨环，所以仅仅靠保留数据点的个数比重和密度等准则是难以保证其加工精度和效率的。点的离散曲率能够较好地反映曲线的真实形状特征，而且计算简单直观，可以很容易地设定曲率阈值，并和加工精度建立联系，所以，本文算法首先引入点云数据的离散曲率作为判断和删除冗余点的依据。

3.1 冗余点分类

根据轮廓环等距和刀位轨迹规划要求，二维截面轮廓中的冗余点主要分为以下两类：

（1）中间冗余点。如果图3中的αi→180°，即称为中间冗余点。

（2）线段较短，短到加工精度以内。

如图6所示，点2、点9与相邻两点组成的线段所成的角接近180°，点4和相邻点在同一条直线上，这3个点统称为第一类冗余点；点5、点11、点12和相邻两点的距离小于加工精度值，为第二类冗余点。

图6 冗余点分类实例

3.2 冗余点判断准则

为满足加工精度要求，减少优化误差，本文算法规定冗余点判断准则包括基于加工精度的粗判断和基于曲率阈值的精判断两步。其判断过程如下：首先把h和加工精度δ的关系作为粗判断条件，即如果hi＞δ，则说明该点是特征点，不能去除，因此不需进行精判断。若hi≤δ，则进入精判断。在精判断中，为了有效地去除图6中的两类冗余点，引入曲率阈值ε的概念，若｜Ki｜＞ε，则说明该点是特征点，不能去除；反之，则说明该点是冗余点，需要去除。可见，使用粗精判断相结合的准则，图6中的两类冗余点就可以合并为一类来处理。

3.3 累积曲率的引入和计算

二维截面轮廓是由一系列小直线段组成的曲线，当点的离散曲率变化比较缓慢时，采用3.2节的方法有可能会去除一系列连续的点，从而误删除有效的特征点，产生过大的累积误差，造成错误的结果。图7为基于离散曲率的点的简化图示。图7a为原始图，其封闭轮廓半环由11个点组成，除了起点和终点外，其他点分布均匀，即K2＝K3＝ …＝K10，假设该曲率值满足｜Ki｜＜ε，如果采用直接判断离散曲率的方法，则会出现图4的简化效果。图7b为去除了前3个点的效果图，由图可见，这时曲线已经偏移了原始轮廓。当遍历到最后一个点时（图7c），原始轮廓就变成了只有起点和终点组成的一条线段，中间所有的点都被删除了，简化结果明显错误。

图7 基于离散曲率的点的简化图示

为了避免这种由于累积误差导致的简化错误，本文提出累积曲率ρ的概念，把离散曲率转化为累积曲率。当某个点确定为冗余点被去除时，其曲率累加到累积曲率中，当计算下一个点时，直接用累积曲率去和曲率阈值作比较，直到该累积曲率值满足点被去除的条件时，ρ清零。在累加曲率过程中，曲率要加入符号位，和曲率阈值比较时，取绝对值对比。

表1是图7a中11个点云数据的离散曲率和累积曲率计算结果，删除标记默认值设为“1”，累积曲率默认值设为“0”，当某点被删除时，置删除标记位为“0”，累积曲率为累积曲率和当前点离散曲率的累加，直到累积曲率值超过曲率阈值时，当前点保留，累积曲率清零。如图8所示，设当前曲率阈值为0.1，若只计算离散曲率，则点2到点9全部被删除，若用累积曲率代替离散曲率，则只有点2、4、6、8、10被删除。

表1 离散曲率和累积曲率的关系

3.4 算法描述

本文提出的基于累积曲率的点云数据简化算法的具体实现步骤如下：

（1）导入二维截面轮廓点云原始数据Pi（i＝1，2，…，n），设累积曲率初始值ρ＝0。

图8 离散曲率与累积曲率

（2）规定起始点为P1，此点定义为定点，定义i的初始值为i＝2，分别计算Δi，hi，li±1，Ki，ρ＝ρ＋Ki。

（3）比较hi和δ的关系，若hi≥δ，βi＝0，ρ＝0，i←i＋1，转（2），否则，转（4）。

（4）比较｜ρ｜和ε的关系，若｜ρ｜＞ε，则该点不是冗余点，保留，βi＝0，ρ＝0，i←i＋1，转（2），否则，转（5）。

（5）去除冗余点Pi，βi＝1。i←i＋1，转（2）。

（6）i＝n时，遍历到最后一个点，程序结束。

4 差分差商自适应光顺算法

简化处理减少了二维截面轮廓点云数据的数量，减少了轮廓曲线在等距、连接、过渡等几何处理中可能出现的一系列不规则问题，有效提高了刀位轨迹生成的速度和质量，降低了实际加工过程中的空走刀、刀具振动等频率。但是，仅仅从数量上对点云数据进行优化是不够的，因为点云数据的分布形态对刀轨生成和实际加工还有重要影响。例如，如果点云数据分布曲率波动较大，很容易造成等距等几何处理后点云排列不均匀，容易产生多处自交等不规则现象，从而不能保证刀轨生成的质量。因此，还需要对简化后的二维截面轮廓点云数据做光顺处理，以使轮廓在加工精度范围内满足光顺要求。

差分差商方法主要是针对点云数据的分布均匀与否来使用的。由于加工对象的形态各异，二维截面轮廓中点云数据的均布程度都是未知的，结合差分差商在曲线光顺中的应用，本文提出差分差商自适应光顺算法。记

式中，Ln为由P1到Pn的累积弦长为P1到Pn的累积弦长的平均值。

对于Pi（i＝0，1，…，n），记

式中，σi为定义的点云数据均布系数，σi越小，说明型值点分布越均匀。

加入判断条件，若σi＜σ，则采用差分法计算点云数据的新位置；若σi≥σ，则采用差商法计算点云数据的新位置。其中，σ为定义的均布因子，本文取σi的中值作为σ。

二维截面轮廓点云数据差分差商自适应光顺算法的流程如图9所示。首先给定点云数据并定义点的坐标Pi，计算相邻点的弦长li、所有点的累积弦长Ln以及平均值，然后计算各点的均布系数σi，并依据此均布系数确定所有点的均布因子σ，比较σi和σ，若σi＜σ，代入式（3）计算新位置点，反之，代入式（4）计算新位置点，并通过式（5）求得差分或差商计算后新位置点的坐标P″i。最后代入式（6），判断新点是否满足误差带容差δ，如果不满足，修改光顺算子α，重新代入式（5），直到满足给定的光顺容差，程序结束。

图9 差分差商自适应算法流程图

5 实例分析

本文算法已应用于Dental Engineer软件中，用图10中的磨牙冠模型（12.2mm×11.4mm×5.7mm）验证该算法的有效性。分别取其3个不同高度的二维截面轮廓为研究对象，对其上的点进行优化，并分别绘制了原始轮廓、优化效果图、10次等距后的效果图，以及简化及优化后点的离散曲率输出图。表2分别记录了3个轮廓环优化前后点的个数、等距一次产生的自交点及等距10次所用的时间。

在图11的二维截面轮廓的简化和等距效果图中，图11a－1、图11b－1、图11c－1分别是磨牙冠模型z＝2.136、z＝2.413、z＝1.986的原始二维截面轮廓图，图11a－2、图11b－2、图11c－2分别是ε＝0.01时的优化效果图，可见，优化后点的数目明显减少，曲线变得更加光滑。图11a－3、图11b－3、图11c－3分别是原始轮廓优化后连续等距10次之后对应的效果图，可见，优化等距后的轮廓质量明显提高，刀轨环变得更加光滑。

图10 磨牙冠模型及3个截平面

图11 二维截面轮廓的简化和等距效果图

图12 二维截面轮廓简化、优化离散曲率输出图

图12为二维截面轮廓简化、优化离散曲率输出图。图12a－1、图12b－1、图12c－1分别是3个二维截面轮廓在ε＝0.02时的优化效果图，通过和图11对比可以看出，ε越小，删除的冗余点越少。图12a－2、图12b－2、图12c－2分别是3个轮廓简化及优化后点的离散曲率输出图，可见：优化算法使点的离散曲率变化更加均匀，多余奇异点和拐点基本消失，曲率极值点减少，保证了轮廓数据二阶导数的连续性，从而确保了加工过程中走刀的稳定性和光顺性。

令α＝0.1，δ＝0.15，ε＝0.02，表2示出了优化前后各轮廓环点的数目、等距1次需要去除的自交点的个数和等距10次的时间。由表2的试验结果可见：通过优化处理，使得轮廓环点云的数目明显减少，需要去除的自交点的数目也减少，节省约一半的等距时间，有效提高了等距等后续处理的效率，从而提高了刀轨环的生成和加工效率。

表2 优化前后二维截面轮廓数据对比

6 结论

本文提出了基于累积曲率的离散刀轨环预处理优化算法。该算法将累积曲率的简化分为粗简化和精简化两步，粗简化保证简化不超出规定的误差带，精简化则确保有效去除冗余点，避免因累积误差导致的非冗余点的误删除。简化算法使二维截面轮廓点的数量大大减少，节省了刀位轨迹生成时间，提高了加工效率。差分差商自适应光顺算法在误差带控制范围内实现二维截面轮廓点的离散曲率均匀化，并保证其二阶导数的连续性，减少了等距过程中自交点分布不均等不规则现象，缩短了计算及处理自交区域的时间，提高了刀位轨迹的质量，改善了产品的加工性能。

［1］Park S C.Sculptured Surface Machining Using Triangular Mesh Slicing［J］.Computer Aided Design，2004，36（3）：279－288.

［2］Reb Y F，Yan H T，Lee Y S.Clean－up Tool Path Generation by Contraction Tool Method for Machining Complex Polyhedral Models［J］.Computers in Industry，2004，54（1）：17－33.

［3］Teh C H，Chin R T.On the Detection of Dominant Points on Digital Curves［J］.IEEE Trans.Pattern Anal.Mach.Intell.，1989，11（8）：859－872.

［4］Hamann B，Chen J.Data Point Selection for Piecewise Linear Curve Approximation［J］.Journal of Computer Aided Geometric Design，1994，11：289－301.

［5］Martin R R，Stroud I A，Marshall A D.Data Reduction for Reverse Engineering RECCAD Deliverable Document 1.COPERUNICUS Project［J］.Computer and Automation Institute of Hungarian Academy of Science，1996，1068：63－69.

［6］刘德平，陈建军.逆向工程中数据精简技术的研究［J］.西安电子科技大学学报（自然科学版），2008，35（2）：334－339.

［7］Liu L，Tai C L，Ji Z，et al.Non－iterative Approach for Globle Meth Optimization［J］.Computer Aided Design，2007，39（9）：772－782.

［8］Hildebrandt K，Polthier K.Constraint－based Fairing of Surface Meshes［C］//Proceedings of the Fifth Eurographics Symposium on Geometry Processing.Barcelona，2007：50－62.

［9］Liu G H，Wong Y S，Zhang Y F.Adaptive Fairing of Digitized Point Data with Discrete Curvature［J］.Computer Aided Design，2002，34：309－320.

［10］Farin G，Spaidis N.Curvature and the Fairness of Curves and Surfaces［J］.IEEE Computer Graphics Applications，1989，9（2）：52－57.

［11］Renz W.Interactive Smoothing of Digitized Point Data［J］.Computer Aided Design，1982，14：267－269.