基于曲率特征的三维模型对称性检测技术
摘 要： 针对三维模型上对称性提取效率不高的问题，在此提出了一种改进的对称性检测方法，该方法通过分析和计算三维模型上曲面的曲率特征，得到模型表面上凸起或者凹陷的褶皱线，再将褶皱线作为分析检测曲面对称性的依据，结合变换空间聚类的方法，计算得到模型整体的对称变换。实验结果表明，该方法通过以简单的褶皱线来描述和刻画模型上的复杂曲面，不但能够得到正确的对称变换，还能够大幅减少算法的计算量，加速提取过程。
关键词： 三维模型处理； 几何建模； 对称性检测； 三维模型曲率特征； 聚类方法
中图分类号： TN964?34 文献标识码： A 文章编号： 1004?373X（2013）08?0089?04
随着CAD及网络技术的快速发展，三维模型的创建和应用越来越广泛，人们每天都在不断地通过手工或3D扫描仪来构建产生虚拟的现实场景及物体，并把它们使用在不同的应用中，例如建筑设计，游戏娱乐，医学分析等。这些三维模型本身大部分都具有对称性的特征，比如建筑，人体，飞机，家具等，这些物体的整体或者部分符合各种各样的对称关系，如直升飞机螺旋桨叶片之间符合45°的旋转关系，汽车车轮之间符合位移关系，大部分的建筑物符合关于中轴对称的反射关系。对称性检测技术就是用来探测这些对称性关系，从而使得人们能够利用这些对称性信息到不同的应用上，例如加速模型编辑过程，模型压缩等。
1 相关工作
近年来，国内外对于对称性检测已经有了不少的研究，大体上分为2种，一种是针对计算模型的内部的（Intrinsic）对称性；另有一种是针对外部的（Extrinsic）对称性。内部的对称性是指模型本身固有的变换特征，并且这种特征不随模型的姿态改变。像Xu等人提出的针对三维模型内部的反射对称变换的探测算法[1]。能够检测出模型上固有的，与模型本身的姿态无关的反射变换，通过在曲面上采样一些点对，使用投票策略来得到一个平滑的标量场，以此来衡量和估算模型上的点是模型固有反射对称轴（Intrinsic Reflectional Symmetry Axis，IRSA）的可能性。然后再通过Grass?fire算法来迭代计算最终的IRSA以及与其对应的曲面。这种方法虽然能够计算与模型姿态无关的对称性，但是仅仅对反射变换有效，且非常耗时。而相对而言，模型的外部对称性的定义就较为严格，它是指模型上近似符合某种变换的曲面特征。变换种类可以是旋转、位移和反射。Mitra等人[2]就曾提出了一种既能抽取反射对称，也能同时抽取旋转和位移对称变换的方法，它通过在变换空间中聚类采样点对，使得显著的变换被抽取，然后再扩散得到对称曲面。同样，Berner和Bokeloh等人[3?4]从另一种角度，使用滑动分析的方法来提取曲面特征，然后提出了一种基于图（Graph）的匹配算法来获得对称性，但是算法整体效率不高。上述的方法虽然各有特点，但是都不高效。为了得到一个较为高效的算法，在此使用曲率作为入手点，分析曲面特征，计算得到曲面上凹陷或凸起的线条位置，称为褶皱线。再利用基于变换空间聚类的对称性检测方法[2]，在褶皱线的基础上分析曲面特征，抽取对称性。这样以“线”代“面”，能极大地减少计算量，提高算法性能，并且能够提升算法对于相对细小的对称关系的敏感性。
2 基于曲率特征的褶皱线提取
曲率是模型表面的基本特征之一，它描述了三维曲面的局部特性，一定程度地刻画和表达了三维曲面，且具有在任何刚性变换下保持不变的特性，因此，在此可以借助这种不变性，来分析和获得相似曲面间的对称变换。通过分析曲率的大小，可以得到一些高曲率或者低曲率的区域，称为特征区域，然后对特征区域进行收缩操作，通过不断收缩，就能得到一些褶皱线，以褶皱线作为对称性抽取的依据能够大大较少计算量，因为大部分的特征不明显的顶点数据已经被过滤掉。2.1节与2.2节将着重介绍特征区域及褶皱线的计算方法。
2.1 特征区域计算
特征区域是模型上凹陷或者隆起的区域，要获得大致的特征区域，首先要估算模型上每个点的曲率大小。对于一个输入的含有N个顶点的三维模型，可以使用文献[5]提出的方法，计算出模型上任意一点Vi的曲率值[Kimax] （可以根据应用背景选取是最大曲率还是最小曲率，或者两者皆可，本文以最大曲率为例）。但是得到的曲率值可能会含有高频的噪声，从而使曲率相对于实际值偏大或者偏小，这时就需要做一些针对噪声的处理。通过大量实验证明，选用均值滤波能够很好的过滤噪声。即对于Vi，其修正后的曲率值[Kimax]等于与Vi相邻的顶点的曲率的平均值。定义集合[F=fi1≤i≤N]表示特征区域，fi与顶点Vi一一对应，且[fi∈0，1]，表示顶点Vi是否被包含在特征区域中。则F中最初的fi可以使用如下公式计算得到：
[fi=1，Kimax>ε0，Kimax≤ε]
式中[ε]为阀值参数，通过实验发现，当[ε=0.8Kmax，][Kmax=max1≤i≤NKimax]时，能得到较好的效果。特征区域的计算结果如图1所示。
2.2 褶皱线提取
直观上讲，褶皱线可以近似地看成模型上的棱，角，凹槽的中心位置，它可以在一定程度上表达曲面的特征。通过前面的方法，已经可以得到一个大致特征区域，接下来可以通过不断地收缩特征区域来最终得到褶皱线。但是，如果仅仅简单的一层层消去特征区域的外层顶点，有可能破坏整个特征区域的连通性，从而使得到的褶皱线不连续，这样的不连续的褶皱线并不能很好地保留曲面的特征。因此，为了保证收缩后特征区域的连通性，对于模型上的每一个顶点Vi，定义：
[Ci=k=0ni-1fuik-fui（k+1）modni]
称Ci为顶点Vi的复杂度，其中ni表示Vi周围一圈点的个数，[uik]表示这一圈点中第k个顶点的编号。对于在特征区域边界上的顶点，其复杂度为2，而在特征区域内部的顶点，由于其邻接的顶点也在特征区域内部，其复杂度为0。当特征区域收缩成线后，线上的顶点的复杂度为4，若该点为若干条线的交点，则其复杂度大于4。注意到如果特征区域收缩成线状的话，其复杂度会变成4，但是如果简单的删掉复杂度等于2的边界点，却会造成线状特征区域在端点处产生退化，因为端点的复杂度也是2。因此，需要定义如下两个辅助集合：
[O=vk=0nvfuvk=nv]
邻接顶点都在特征区域上的点都包含在集合O中，如果特征区域收缩成线状的话，线上必定没有集合O中的点。通过集合O可以得到其中每一个点的邻接顶点的并集：
[D=v?p∈O，且v∈up]
因此，线状的特征区域的端点必然不包含在D中。所以，在收缩的时候，只需不断地从特征区域上删除在D中且复杂度等于2的顶点，即可维持特征区域的连通性。因此计算褶皱线的算法步骤为：
（1）计算特征区域上顶点的复杂度以及集合D。
（2）删除集合D中复杂度为2的边界顶点，若没有顶点被删除，则完成计算，否则跳至步骤（1）。
通过上述步骤，即可得到如图2所示的褶皱线，图中黄色的线条即为图1所示的特征区域通过上述方法得到的褶皱线。
3 对称性检测
褶皱线反应了曲面的一些固有的特征，可以利用这种特征来分析曲面间的对称关系，结合文献[2]的对称性检测方法，通过分析每两个褶皱线上的点之间的变换，使用类似投票的方法，让每个点对在变换空间中对其所反映的变换投票，最后票数最多的变换即为模型上显著的对称变换。
3.1 配对及变换的计算
与文献[2]不同，本文使用褶皱线上所有的点作为采样点，这些点构成了点集P。对于P中任意一个点Pi，它的2个主曲率，即最大曲率值和最小曲率值分别为[ki1]和[ki2]，它们反应了Pi点的局部几何特征，定义其二维特征向量[Ki=ki1，ki2]，当Pi与[P]中另外一个点Pj（[i≠j]）进行配对时，[Ki]与[Kj]在二维平面空间中的欧氏距离是否在限定的范围之内，如若符合要求，则计算Pi到Pj的变换矩阵。使用文献[6]所述的近似最近邻方法，可以在O（|P|log |P|）的时间内完成上述步骤。任意两点Pi和Pj的变换矩阵是一个4×4的矩阵，它可以由旋转矩阵rij和位移向量tij合成得到，因而Pi和Pj之间的变换tij可以表示成一个六元组[Tij=rxij，ryij，rzij，txij，tyij，tzij]，其中[rxij，ryij，rzij]表示旋转矩阵rij的3个欧拉角，[txij，tyij，tzij]是位移向量tij的3个分量。使用Pi的主曲率方向[ci1]，[ci2]和法向[ni=ci1×ci2]，构造Pi上的局部坐标系[ci1，ci2，ni]，同样构造Pj上的局部坐标系[cj1，cj2，nj]，这两个局部坐标系间的变换即为旋转矩阵rij，然后可计算[tij=Pj-rijPi]。对于反射变换的计算，可以先把Pi反射到以固定平面为轴的镜面空间中，再计算其变换，本文以x=0为固定轴进行反射。这样，当每一对Pi和Pj之间的变换Tij都计算出来后，即可得到一个包含了[P2]个点的6维变换空间。
3.2 聚类
对于每一个Tij，都是由两个局部几何特征相似的点计算得到的，且Tij是?中的点。如果模型上两块曲面上的所有顶点都完美的符合同一个变换，则这些顶点所构成的变换肯定都落在?中的同一个点上。但是实际上，由于曲率估算的误差，噪声的影响等因素的存在，这些点对所反映的变换不可能是完全相同的。因此需要使用聚类的方法，来辨别出比较显著的变换的。由于分类后的类的个数是未知的，所以无法使用需要先验知识的分类方法，如K?means[7]。因此，采用了文献[8]中所提出的mean?shift方法，并使用对称的Epanechnikov核。该方法把空间?中的点分成了若干个类，对于任意一个类Gk，其对应的变换为Tk，即该类的局部密度最大值在空间?中的坐标。
3.3 对称曲面的验证
通过聚类，可以得到由若干个点对组成的Gk，且这些点对都近似符合变换Tk。利用这些点对，可以使用双向扩散的方法，来计算出对称的曲面。对于Gk中任意一个点对（Pi，Pj），以Pj为参照，保留Pi周围的一圈顶点中，经过Tk变换后，与Pj的距离小于某一给定阈值的顶点，记这些被保留下来的顶点集合为Si。同样的，对于Pj，以Si为参照，判断其周围一圈顶点，产生集合Sj。不断地扩大Si和Sj，直至涵盖了用户将要操作的感兴趣区域为止。上述方法涉及到空间中最近距离的查询问题，可以维护一个KD树，这样每次查询和新增顶点的操作都可以在O（log |S|）时间内解决。由于聚类出来的变换Tk是近似的，不准确的，所以还需要逐步的求精，在不断扩散得到Si和Sj的过程中，可以借助迭代最近点对算法[9]（Iterative Closest Points Algorithm）不断的迭代求精Tk。
4 实验结果
为了验证本文改进方法的可行性以及有效性，从算法的运行效果以及运行效率2个方面进行了实验，并与文献[2]中的方法进行了对比。实验数据采用了PSB标准模型库[10]以及互联网上的一些模型，分别针对位移变换，旋转变换，反射变换3种对称性关系做了提取，得到了模型上的所有对称变换，其中最显著的对称变换所对应的对称曲面如图3所示。
第1列为输入的模型，第2列中红色区域和黄色区域表示对称表面。其中上，中，下3行分别对应了旋转变换，位移变换和反射变换。本文所有试验均在一台具有4 GB内存，Intel Core I7 处理器的便携式计算机上进行。运行环境是Ubuntu 10.10。同时为了验证本文提出的方法的高效性，实现了一种基于变换空间聚类（Transform Space Clustering，TSC）且未使用任何曲面特征的对称性检测算法[2]，简称TSC算法，并将该算法与本文中的改进算法进行了对比。实验记录了本文算法和TSC算法得到模型上全部对称变换所用的时间，以及本文算法提取褶皱线的耗时，另外，为了验证本文所使用的褶皱线特征的有效性，统计了TSC算法为了得到正确结果所需要的最小数目的采样点，以及本文所使用的采样点数目，即所有褶皱线上的顶点数目，实验结果如表1所示。
从表1可以看出，本文算法相比TSC算法，虽然增加了褶皱线提取的耗时，但是褶皱线能够大幅减少了采样点数目，使得后续的配对和聚类过程的耗时大幅降低，从而使得算法的整体效率有明显提升。
5 结 语
本文提出了一种基于曲率特征的对称变换的提取算法，它使用褶皱线来描述三维模型上的复杂曲面，大大减少了对称变换提取时所需要处理的顶点数，从而提升了算法效率。它相对于其他的对称性检测方法，具有更广的适用性，更能够给用户带来高效的体验。
参考文献
[1] XU Kai， ZHANG Hao， TAGLIASACCHI Andrea， et al. Partial intrinsic reflectional symmetry of 3D shapes[C]// ACM SIGGRAPH Asia 2009 papers， Yokohama， Japan： CM SIGGRAPH， 2009： 16?19.
[2] MITRA Niloy J， GUIBAS Leonidas J， PAULY Mark. Partial and approximate symmetry detection for 3D Geometry [J]. ACM Transactions on Graphics， 2006 25 （3）： 30?36.
[3] BOKELOH M， BERNER A， WAND M， et al. Symmetry detection using feature lines [C]// Comput. Graph. Forum. [S.l.]： CCF， 2009： 697?706.
[4] BERNER A， BOKELOH M， WAND M， et al. A graph?based approach to symmetry detection [J]. Volume Graphics， 2008， 40： 1?8.
[5] COHEN?STEINER D， MORVAN J. Restricted delaunay triangulations and normal cycle [C]// Symposium on Computational Geometry. [S.l.]： SCG， 2003： 312?321.
[6] MOUNT D M. ANN Programming Manual [M]. Maryland： Department of Computer Science， University of Maryland， College Park， 1998.
[7] HARTIGAN J A， WONG M A. A k?means clustering algorithm， applied statistics [M]. [S.l.]： [s.n.]， 1979.
[8] COMANICIU D， MEER P. Mean shift： a robust approach toward feature space analysis [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence， 2002， 24（5）： 603?619.
[9] RUSINKIEWICZ S， LEVOY M. Efficient variants of the ICP algorithm [C]// Proceedings of 2001 Third International Conference on 3?D Digital Imaging and Modeling. Quebec City， Canada： [s.n.]， 2001： 145?152.
[10] SHILANE Philip， MIN Patrick， KAZHDAN Michael， el al. The Princeton shape benchmark [C]// Proceedings of 2004 Digital Object Identifier. Genova： DOI， 2004： 167?178.
