基于改进圆球覆盖算法的阶梯状伪影去除方法

陈 中，段黎明，吴志芳

（1.重庆大学 光电技术及系统教育部重点实验室ICT 研究中心，重庆 400030；2.重庆大学 机械工程学院，重庆 400030；3.机械工业第三设计研究院，重庆 400039）

0 引言

在逆向工程领域，工业CT（Industrial Computed Tomography，ICT）技术是一种先进的非接触测量方法，通过该方法可以得到工件的外表面和复杂内腔的测量数据。由于受工业CT 扫描成本和扫描工艺的限制，工件的扫描间距较大，导致断层切片图像的层间距离远大于切片图像的平面分辨率，使用这些CT 断层序列图像重建工件的三维表面网格模型时，会产生明显的阶梯状伪影，影响了后续的特征提取及三维实体模型的逆向重构［1-2］。

目前，消除阶梯状伪影的主要途径是使用光顺算法对重建的网格模型进行表面平滑处理。然而，光顺算法一般会产生过光顺现象，例如拉普拉斯［3］（Laplacian）光顺算法不能保持模型的精度，也不能有效去除阶梯状伪影。保特征的光顺算法［4-6］在去除表面模型噪声的同时能有效保持模型特征，但是该算法会将阶梯状伪影作为模型特征加以保留。文献［7］针对医学表面模型提出一种基于微分域的自适应光顺方法，该方法可以去除噪声和表面伪影，但是效果不明显。以上光顺算法都是在网格模型生成后对其进行网格平滑处理，而且都不能有效去除阶梯状伪影。近些年，隐式曲面重构光顺网格模型方法的应用日益广泛，这种方法通过局部拟合的方式消除一部分表面网格伪影。文献［8］采用正交多项式作为拟合函数，对多层次单元分解（Multi-level Partition of Unity，MPU）方法进行改进，直接从点云数据中重构出光顺的三维网格模型，但是该方法不能生成封闭的高质量网格模型。文献［9］提出一种点云网格化的集成方法，集成了网格重建和孔洞修补等功能，该方法主要通过构造二次误差函数来拟合点云，从点云直接重建成封闭的高质量网格模型，但是该方法对夹杂在点云中的噪声比较敏感。如果点云中夹杂的噪声较多，则会生成低质量网格模型，甚至导致网格重建失败。文献［10］对文献［9］中的相关方法进行了改进，在结构噪声较少的情况下，通过自适应采样处理分布非常不均匀的点云数据，使得点云法向量方向一致化，从而正确重构曲面。

针对上述问题，本文根据工业CT 逆向数据的特点，提出一种新的去除表面阶梯状伪影的方法，从工业CT 点云数据中直接重建光顺的网格模型。首先采用高斯加权的协方差矩阵估算出工业CT 点云的法矢量，然后使用双边滤波去除点云中的噪声，最后通过二次误差函数拟合点云，使用自适应圆球覆盖方法对点云进行网格化处理。通过实例验证表明，本方法可以从工业CT 点云模型直接生成光顺的封闭网格模型，而且在生成高质量网格模型的同时有效去除了阶梯状伪影，避免了后期的网格优化处理。

1 点云法矢量估算

点云去噪和网格重建需要估算点云的法矢量。目前对点云法矢量估算的方法较多［11-14］，主成分分析（Principal Component Analysis，PCA）法［15］是其中比较经典的方法，但是传统PCA 法的鲁棒性较差，当点云含有较多噪声时，会使法矢量产生偏差，不能正确估算点云法矢量，甚至导致曲面重建失败。产生上述问题的主要原因是在拟合最小二乘平面时，噪声点与采样点使用了相同的权值。CT 图像的噪声一般近似为符合高斯分布的加性噪声［16-17］。在CT 序列图像转化为三维点云的过程中，序列图像含有的高斯噪声也转化为相应的噪声点存在于点云中，点云局部呈现高斯分布。本文采用高斯加权的协方差矩阵估算点云法矢量，可以根据点云中的局部不同点对整体影响程度的不同赋予相应的权值，另外高斯函数对去除服从正态分布的噪声有很好的效果，因此权重函数选择高斯函数，从而有效降低高斯噪声的干扰，比较准确地估算出点云法矢量。

为了准确重建扫描工件的三维表面网格模型，首先建立工业CT 点云的k-邻域，用pj表示为k-邻域中的点表示为k-邻域的质心。本文采用kd-tree作为邻近点查询的数据结构，选用k-邻域查询方式处理非规则模型，其自适应性较强。k 值反映模型的局部性质，考虑到工业CT 点云分布的情况和查询效率，本文选择k=10。

对以下高斯加权的协方差矩阵C 作特征值分解：

式中：ψj为高斯函数，μ为点pi到其k-邻域中所有点的平均距离，dj为点pi到其k-邻域中任一点pj的距离。求解出的三个特征值记为λ3≥λ2≥λ1，该协方差矩阵的最小特征值λ1对应的单位特征向量可以作为点pi的单位法矢量ni。

2 点云去噪

常用去噪技术是通过滤波来实现的［18］，对含有噪声的点云进行滤波处理，可以有效降低点云的随机误差，提高网格模型重建的精度。双边滤波方法最早被用于图像的滤波去噪处理，Fleishman［4］和Jones［5］等将双边滤波法推广到三维网格模型的光顺操作，该方法可以有效去除网格模型的噪声并保持模型的特征。双边滤波法只与滤波点的邻域有关，也适用于点云模型的滤波去噪。通过双边滤波函数计算出点的移动量d，结合点pi的单位法矢量ni，可以得到滤波后的点。本文通过双边滤波方法进行多次迭代能够有效去除工业CT 点云数据中大量的噪声。

式中：k 为距离点pi最近的邻域点数，Wc（x）为空间权函数，Ws（x）为特征保持权函数，σc为切平面的高斯滤波系数，σs为法向的高斯滤波系数。

3 点云网格化

对于切片间距过大引起的阶梯状伪影表面，可以根据表面阶梯情况选择合适的参数，通过自适应圆球将所有层间阶梯进行覆盖，使用局部拟合的方式逼近光顺的表面。工业CT 点云数据是均匀有序分布的［19］，没有非常不均匀的点数据，不会出现覆盖球体集完全分离的情况［20］。对工业CT 点云去噪后，得到新的点云是滤波后点的k-邻域中的任一点。

定义二次误差函数

式中：c为圆球的球心｛c1，…，cm｝∈P′，r为圆球的半径｛r1，…，rm｝，x为辅助点，wi为点云的权值，Gσ（ρ）为类高斯函数，Tq为用户控制圆球的参数。

求解二次误差函数的最小值，可以转换为求解xmin。

式中：E（r）为误差函数，L 为点云P′的包围盒主对角线长度。

通过用户自定义的精度Terr可以计算出xmin。根据‖x-c‖≤r判断，如果xmin在圆球内部，则将xmin作为辅助点，如果xmin在圆球外部，则球心c为辅助点，最后将m 个相交圆球的辅助点连接起来得到三角网格模型。本方法总体流程示意图如图1所示，算法流程如下：

步骤1 输入工业CT 点云数据，用户设定相关参数。

步骤3 点云P 中任意一点pi经过双边滤波后变为pi′，形成新的点云。

步骤4 通过用户自定义的精度，使用二次误差函数拟合去噪后的点云，根据‖x-c‖≤r判断，将m 个相交圆球的辅助点相互连接得到三角网格模型。

4 应用实例

本文在VC++6.0的环境下实现了上述算法，在Intel（R）Core（TM）i3CPU 2.53GHz，2GB内存的计算机上进行了实验。本文算法的时间复杂度为O（nlog n），空间复杂度为O（n），n为输入点云的点数。以汽车轮毂和茶壶的CT 数据重建网格模型为例，介绍使用自适应表面网格重建法去除阶梯状伪影的整个过程。轮毂通过工业CT 扫描获得67张断层序列图像，如图2a所示。轮毂CT 数据的体素大小为0.315mm×0.315mm×3mm，即切片平面xy 像素大小为0.315 mm×0.315 mm，切片层间距z 为3mm。茶壶层间距为3mm，共有59张断层序列图像，如图2b所示。目前针对CT 数据的可视化主要是将切片序列图像直接生成三维模型，一般使用体绘制和面绘制方法。体绘制方法可以比较真实地反映工件的内外结构，面绘制方法可以重建工件的三维表面网格模型。图3所示为轮毂模型和茶壶模型的体绘制效果图，图4所示为轮毂模型和茶壶模型的面绘制效果图。从图3和图4中可以看出，轮毂曲面和茶壶曲面呈阶梯状，但是真实轮毂和茶壶的曲面是光滑的，这是因为扫描间距过大而在重建时产生表面伪影。

本文首先将工业CT 的序列图像转化为三维点云并估算出点云的法矢量，然后通过双边滤波去除点云噪声，接着使用自适应圆球覆盖方法对点云进行网格化处理。自适应调整圆球的大小，曲率变化较大的表面区域圆球半径较小，曲率变化较小的区域圆球半径较大，相交圆球的辅助点相互连接最终形成三角网格模型。

为验证本文方法去除阶梯状表面伪影的效果，与常用的去除阶梯状伪影的方法——拉普拉斯方法和Fleishman网格光顺算法进行对比，待光顺的模型是由工业CT 序列图像通过面绘制MC（marching cubes）算法生成的网格模型。不同方法的对比实验都是以去除阶梯状伪影为准则，伪影去除为上述方法处理的最佳结果。使用拉普拉斯光顺方法对轮毂网格模型经过15次迭代后能去除部分阶梯状伪影，23次迭代后能去除所有阶梯状伪影，但是模型变形非常严重。依次使用15～23之间的次数进行迭代实验，在不破坏原有模型结构的情况下，迭代19次时能取得最佳的去伪影效果；茶壶模型迭代21次时能取得最佳效果。Fleishman方法具有很好的特征保持性，会将阶梯状作为模型特征保留下来，该方法经过2 次迭代，就能够获得最佳的去伪影效果。本文算法参数取值经过验证取值为：去噪迭代次数n=3，σc=3，σs=2，Terr=1.0×10-5，Tq=2。图5所示为不同方法处理轮毂模型后的效果图，图6所示为圆球覆盖算法改进前后轮毂模型的对比效果图，图7所示为不同方法处理茶壶模型后的效果图，图8所示为圆球覆盖算法改进前后茶壶模型的对比效果图。从图5和图7可以看出，拉普拉斯光顺方法虽然经过多次迭代去除了表面阶梯状伪影，但是产生了较为严重的过光顺现象，轮毂模型变形较大且螺孔边缘位置收缩明显，茶壶盖顶部的把手变形严重。Fleishman 方法有较好的特征保持性，轮毂螺孔边缘没有明显变化，茶壶盖顶部的把手也没有明显变形，但是阶梯状伪影被作为特征保留下来。本文方法有效去除了轮毂和茶壶表面的阶梯状伪影，轮毂的螺孔边缘保持较好，茶壶盖顶部的把手没有产生严重变形，但由于受圆球半径选取大小等情况的影响，轮毂边缘有些粗糙。

文献［9］是一种集成的方法。本文方法针对工业CT 点云数据的特点改进了点云法矢量估算方法，并且将双边滤波集成到该方法中，为了直观显示改进后的效果，与文献［9］的方法进行了详细对比。文献［9］方法的参数设置为k=10，Terr=1.0×10-5，Tq=2。图6a～图6c以及图8a～图8c为文献［9］方法的处理结果，图6d～图6f以及图8d～图8f为本文方法的处理结果。虽然两种方法设置的参数k，Terr和Tq相同，输入的点云相同，但是由于受噪声影响，文献［9］方法的圆球覆盖数量较多，最终重建的网格模型表面比较粗糙，网格数量较多，具体对比数据如表1所示。表2所示为四种方法处理后的网格模型与离散化的原始计算机辅助设计（Computer Aided Design，CAD）模型的误差比较，从表2可以看出，采用本文方法重建的模型误差较小。

表1 算法改进前后的数据统计

续表1

表2 各种方法的误差比较 mm

5 结束语

本文实现了一种去除阶梯状伪影的网格重建方法。该方法通过高斯加权的协方差矩阵估算出工业CT 点云的法矢量，引入双边滤波去除点云噪声，使用二次误差函数拟合点云，自适应圆球覆盖对点云进行网格化处理，最终形成三角网格模型。将该方法应用于工程实际，实例表明，本文方法在生成高质量三角网格模型的同时能有效去除阶梯状伪影，并且网格模型精度保持较好。本文方法受圆球半径选取大小等情况的影响，重建的模型边缘有些粗糙，拟集成一种边缘特征增强方法进行改进，这是下一步的研究方向。

［1］HE Qiang，ZHANG Shusheng，BAI Xiaoliang.Robust algorithm for symmetry plane extraction in reverse engineering［J］.Computer Integrated Manufacturing Systems，2012，18（6）：1121-1126（in Chinese）.［贺 强，张树生，白晓亮.逆向工程中一种鲁棒的对称面提取方法［J］.计算机集成制造系统，2012，18（6）：1121-1126.］

［2］DUAN Liming，LIU Yuanbao，WU Zhifang，et al.Method of reconstructing 3-D CAD model based on industrial computed tomography［J］.Computer Integrated Manufacturing Systems，2009，15（3）：479-486（in Chinese）.［段黎明，刘元宝，吴志芳，等.基于工业计算机断层成像技术的三维CAD 模型重构方法［J］.计算机集成制造系统，2009，15（3）：479-486.］

［3］VOLLMER J，MENCL R，MÜLLER H.Improved laplacian smoothing of noisy surface meshes［J］.Computer Graphics Forum，1999，18（3）：131-138.

［4］FLEISHMAN S，DRORI I，COHEN-OR D.Bilateral meshes denoising［C］／／Proceeding of SIGGRAPH 2003.New York，N.Y.，USA：ACM，2003：950-953.

［5］JONES T R，DURAND F，DESBRUN M.Non-iterative feature-preserving mesh smoothing ［C］／／Proceedings of SIGGRAPH 2003.New York，N.Y.，USA：ACM，2003：943-949.

［6］ZHANG Xin，WANG Zhangye，FAN Hanqi，et al.A feature preserving denoising approach for scanned models based on trilateral filtering［J］.Journal of Computer-Aided Design &Computer Graphics，2009，21（7）：936-942（in Chinese）.［张鑫，王章野，范涵奇，等.保特征的三维模型的三边滤波去噪算法［J］.计算机辅助设计与图形学学报，2009，21（7）：936-942.］

［7］WU Jianhuang，LI Yonghong，MA Xin，et al.Medical surface smoothing via adaptive diffusion of differential fields［C］／／Proceedings of the 2011International Conference on Complex Medical Engineering.Washington，D.C.，USA：IEEE，2011：323-327.

［8］WANG Yanan，YANG Bo，ZHAO Xiuyang.Reconstruction of cement particle by contour extraction and adaptive partitioning of unity［J］.Journal of Computer-Aided Design &Computer Graphics，2010，22（1）：152-157（in Chinese）.［王亚楠，杨 波，赵秀阳.基于轮廓提取和自适应单元分解的水泥颗粒三维重构［J］.计算机辅助设计与图形学学报，2010，22（1）：152-157.］

［9］OHTAKE Y，BELYAEV A，SEIDEL H P.An integrating approach to meshing scattered point data［C］／／Proceedings of ACM Symposium on Solid and Physical Modeling.New York，N.Y.，USA：ACM，2005：61-69.

［10］LIU Shengjun，HAN Xuli.Sampling and orienting of scattered points based on adaptive spherical covering［J］.Chinese Journal of Computers，2011，34（3）：489-498（in Chinese）.［刘圣军，韩旭里.基于自适应球覆盖的散乱点采样与定向［J］.计算机学报，2011，34（3）：489-498.］

［11］GUENNEBAUD G，GROSS M.Algebraic point set surfaces［J］.ACM Transactions on Graphics，2007，26（3）：231-239.

［12］YOON M，LEE Y，LEE S，et al.Surface and normal ensembles for surface reconstruction［J］.Computer-Aided Design，2007，39（5）：408-420.

［13］HUANG Hui，LI Dan，ZHANG Hao，et al.Consolidation of unorganized point clouds for surface reconstruction［J］.ACM Transactions on Graphics，2009，28（5）：176.

［14］LI Bao，SCHNABEL R，KLEIN R，et al.Robust normal estimation for point clouds with sharp features［J］.Computers&Graphics，2010，34（2）：94-106.

［15］PAULY M，GROSS M，KOBBELT L P.Efficient simplification of point-sampled surfaces［C］／／Proceedings of IEEE Viualization 2002.Washington，D.C.，USA：IEEE，2002：163-170.

［16］WANG Hao，KANG Xiaodong，LIU Lingling，et al.Medical CT image denoising method based on sampled Contourlet transform［J］.Computer Engineering and Applications，2012，48（27）：150-154（in Chinese）.［王 昊，康晓东，刘玲玲，等.基于非下采样Contourlet变换的医学CT 图像去噪［J］.计算机工程与应用，2012，48（27）：150-154.］

［17］ZHA Fanglong，ZHANG Dinghua，HUANG Kuidong，et al.Subvoxel-accuracy surface detection method based on contour pre-segmentation for computed tomography images［J］.Chinese Journal of Scientific Instrument，2012，33（6）：1308-1314（in Chinese）.［查方龙，张定华，黄魁东，等.基于预分割轮廓的CT 图像亚体素表面检测方法［J］.仪器仪表学报，2012，33（6）：1308-1314.］

［18］LI Jinjiang，ZHANG Caiming，FAN Hui.Point cloud denoising algorithm based on swarm intelligent［J］.Computer Integrated Manufacturing Systems，2011，17（5）：935-945（in Chinese）.［李晋江，张彩明，范 辉.群体智能点云光顺去噪算法［J］.计算机集成制造系统，2011，17（5）：935-945.］

［19］LI Xingdong，LAI Xinghua，LIU Wenli，et al.Extracting point cloud data for jacket of vehicle crash test dummy based on CT images［J］.Chinese Journal of Stereology and Image Analysis，2011，16（3）：254-261（in Chinese）.［李兴东，赖兴华，刘文丽，等.基于CT 图像的汽车碰撞假人上衣模型点云数据提取研究［J］.中国体视学与图像分析，2011，16（3）：254-261.］

［20］LIU Shengjun，WANG C L.Orienting unorganized points for surface reconstruction［J］.Computers &Graphics，2010，34（3）：209-218.