0 引言

近年来，散乱点云的曲面重建技术在逆向工程、计算机辅助设计与制造、医学图像重建、虚拟现实和科学计算可视化等领域得到了广泛的研究与应用^[1-2]。目前，曲面重建方法主要分为三种^[3]：基于Delaunay三角化的曲面重建、隐式曲面重建和区域生长法曲面重建。其中，基于Delaunay三角化的曲面重建方法是唯一具有严格理论证明的网格曲面重建方法，在曲面重建算法中受到广泛关注和研究，其最经典的方法是α-shape算法^[4]、Crust算法^[5]和Cocone算法^[6]。

曲面重建算法依赖于点云类型，不同类型的点云需要不同的曲面重建方法^[7]。通过三维激光扫描仪等设备获取的点云通常是非均匀的散乱点集且仅仅包含点的空间三维坐标信息，这对曲面重建算法提出了苛刻的要求。α-shape算法思路简单，易于实现，但要求采样密度均匀，当采样点分布不均匀时，很难选择合理的α值重建出物体表面。Xu等^[8]和Teichmann等^[9]依据局部点云密度自适应地调整α值,但点云密度估计质量没有评估标准且多元自适应过程有待进一步研究。Giesen等^[10]根据各个单纯形的α值序列α_pⁱ和采样点的局部特征尺寸，构造出内部α尺度参数_pⁱ，实现自适应的α-shape重建，但重建结果存在大量孔洞且重建效率降低。Jiang等^[11]基于滚球法(Ball Pivoting Algorithm,BPA)，根据采样点间距对服装表面网格进行重建，极大地减少了α-shape算法重建结果中的不期望三角面片，但没有实现自适应重建。Lafarge等^[12]给出修正版本的alphashape3d库，不仅实现了α-shape算法，还优化了三维图形显示和重要特征计算的功能，但没有考虑如何处理非均匀点集。胡丝兰等^[13]使用可变半径改进BPA算法，使得算法可以处理不均匀点云数据，但是可变半径值是重建前手动设置的。孙殿柱等^[2]用增益优化后的局部样本自适应地调整α值，可显著减少重建结果中的孔洞且算法效率与主流Delaunay算法相当，但不适用于海量点云的曲面重建。上述非均匀点集的α-shape曲面重建算法都存在一些局限性，有必要对α-shape算法作进一步的研究。

针对α-shape算法不适用于非均匀点集曲面重建的问题，本文在文献^{[4, 10]}的基础上，提出了一种基于局部特征尺寸的散乱点云的自适应α-shape算法。该算法可以降低噪声对局部特征尺寸估计的影响，简化后的点云减少了曲面重建的耗时，并且可以实现非均匀点集的自适应曲面重建。

1 基本定义

局部特征尺寸计算是本文算法的基础，下面给出与局部特征尺寸密切相关的几个定义^{[3, 5]}：

Voronoi图：三维空间中的点集P={p_i|i=1,2,…,n},距离点p_i最近的点集V_i={q||qp_i|≤|qp_j|,j≠i,p_i,p_j∈P}，其中，|qp|表示点q到点p的欧氏距离。称V_i为点p_i的Voronoi单元，点集P中的所有点的Voronoi单元的并集为P的Voronoi图。

曲面中轴(Medial Axis,MA)：距离空间曲面S至少有两个最近点的点集称为S的曲面中轴MA。

极点：采样点p对应的Voronoi单元的两个最远的且位于空间曲面S两侧的顶点称为p的极点。

局部特征尺寸(Local Feature Size,LFS)：空间曲面S上的任意点p到曲面中轴的最近距离f(p)称为曲面S在采样点p处的局部特征尺寸。

ε采样：曲面S的采样点集为P，对于S上任意点p，P中都存在点q满足‖p－q‖≤εf(p)，则称点集P满足ε采样条件。

2 自适应α-shape重建算法

自适应的α-shape重建算法根据局部区域采样点的细节信息自动调节参数α的大小。本文先计算采样点的局部特征尺寸和相邻采样点间距，再依据局部特征尺寸简化点云，最后自适应重建曲面。

2.1 α-shape算法

1994年，Edelsbrunner等^[4]第一次给出了α-shape的通用概念，并给出点集α-shape族的一种计算算法。α-shape算法先对散乱点集进行Delaunay三角剖分，然后对剖分结果中的每个单纯形(四面体、三角面、边、顶点)分别计算其属于α-shape的取值区间，若α值位于该取值区间内则保留该单纯形，若α值不在该取值区间内则删除该单纯形。

如图 1所示，α-shape算法的参数α不能自适应地确定，重建结果随着α值的改变而改变，α值取值偏大会保留过多的三角形、四面体和边，α值取值偏小会出现孔洞。要想获取能够重建出物体表面的合适的α值需要经过多次尝试，但是对于采样密度不均匀的采样点，只取一个固定的α值不能重建出高质量的曲面。

2.2 局部特征尺寸估计

在曲面重建过程中，因曲面是未知的，故严格意义上的曲面中轴计算十分困难，同样直接计算采样点的局部特征尺寸也非常困难。当点云不存在噪声时，Amenta等^[5]发现利用采样点的极点可以较好地逼近曲面中轴。但是，单点计算的极点容易受点云噪声影响。当点云存在噪声时，部分极点会偏离曲面中轴而靠近曲面，以致无法准确估计曲面中轴，如图 2(a)。本文用k-邻近点为整体得到的公共极点逼近曲面中轴，该方法在点云存在噪声时依然能准确地估计出曲面中轴^[14]，如图 2(b)所示。

本文计算采样点的局部特征尺寸的流程为：

1) 对点集P进行Delaunay三角化；

2) 对于P中任一点p，计算出其k-邻近点N_p；

3) 计算正极点，即N_p中各点对应的Voronoi单元的最远的且在曲面S一侧的顶点v_p1；

4) 计算负极点，即N_p中各点对应的Voronoi单元的最远的且在曲面S另一侧的点v_p2。

5) 提取出所有的负极点组成集合Q。

集合Q就是曲面中轴的估计，采样点到Q的最近距离为曲面在该点处的近似局部特征尺寸。

2.3 点云简化

三维激光扫描仪获取的点云数据十分庞大，这些数据不仅占用计算机大量的物理内存而且会增加曲面重建耗时，故需要对点云数据进行简化。点云数据简化的关键是在简化点云的同时，能够最大限度地保留点云数据的细节特征^[15]，即在细节特征信息丰富区域删除少量的点；在细节特征信息匮乏区域删除较多的点。

采样点p的局部特征尺寸f(p)反映了曲面在该点处的细节特征变化程度，在特征变化越显著的区域f(p)越小，在特征变化不明显的区域f(p)越大。曲面上与采样点p距离小于εf(p)的区域的细节特征可通过该点充分地表达^[14]，可删除该区域内的其余内部点，即可以根据曲面的局部特征尺寸对点云进行降采样。

2.4 自适应曲面重建

对所有简化后的点云重新进行Delaunay三角剖分，对剖分结果中的每一个三角面片进行判断，设三角面片的外接球半径为r，若r≤α,则保留该三角面片；否则，删除该三角面片。

2.4.1 α值的自适应调整

如图 3所示，Delaunay三角剖分结果中，物体表面的每一个三角面片pqr可近似为等边三角形，且在该三角面片的微小区域内，所有点的局部特征尺寸近似相等，因此可以用三角面片pqr三个顶点的局部特征尺寸的均值作为该三角面片的局部特征尺寸f(t)的估值，并用该均值作为三角面片三个顶点的新的近似局部特征尺寸，即f(p)=f(q)=f(r)=f(t)。本文根据每一个Delaunay三角面片的局部特征尺寸自适应调整α值的大小。

由几何知识易知，△pqr的最小外接球半径等于其最小外接圆半径。如图 3所示，正△pqr的最小外接球半径为：

$r={l}/{\sqrt{3}}\;$

(1)

式中，l为Δpqr的边长，亦即简化后相邻采样点的间距。如图 4所示，对原始点云数据进行简化时，曲面上与相邻采样点p、q距离不大于εf(t)的采样点都被删除，简化后相邻采样点间距的取值范围为:

$\varepsilon f(t)＜l＜2\varepsilon f(t)+d$

(2)

其中d为原始点云相邻两个采样点之间的间距。

进一步分析可知：在曲面平滑处，三角面片pqr的局部特征尺寸相对较大，d远小于εf(t)。在曲面高曲率处，三角面片pqr的局部特征尺寸会很小，甚至d远大于εf(t)，此时l远大于2εf(t)，以致可以忽略εf(t)的影响。因此在局部特征尺寸相对较大的曲面处可以忽略d的影响，但是在局部特征尺寸很小的曲面处若忽略d的影响会导致最后计算出的结果出现大量的孔洞，此时d对l的影响必须考虑。本文用采样点p到其k-邻近点内的其他点的距离的均值作为采样点p的相邻采样点间距d_p，在△pqr所在的区域内的原始点云相邻点间距大致相等，因此用△pqr的三个顶点的采样点间距的均值作为式中的d值是合理的。

由式(1)~(2)联立可知:

${\varepsilon f(t)}/{\sqrt{3}}\;＜r＜{\left( 2\varepsilon f(t)+d \right)}/{\sqrt{3}}\;$

(3)

因为当r≤α,则保留该三角面片，为尽可能减少重建输出的α-shape表面的孔洞，α值取上界，即:

$\alpha ={\left( 2\varepsilon f(t)+d \right)}/{\sqrt{3}}\;$

(4)

2.4.2 三角面片的外接球半径r的确定

对于每个Delaunay三角面片，其顶点坐标已知，根据简单的几何知识很容易求出其最小外接圆半径r。具体步骤为：

1) 计算三角面片pqr的三边边长a、b、c和周长C；

2) 利用海伦公式$S=\sqrt{u\left( u-a \right)\left( u-b \right)\left( u-c \right)}$计算三角面片的面积S，其中u为周长C的一半；

3) 利用公式r=abc/(4S)计算三角面片的最小外接球半径r。

3 实验结果及讨论

为测试改进算法效果，选择下载于http://www.cc.gatech.edu/projects/large_models/网站的点云数据bunny、blade和dragon为实验数据，在Windows 7系统Intel Core2 i5-2450M CPU 2.50 GHz硬件环境下，结合计算几何算法库(Computational Geometry Algorithms Library,CGAL)^[16]实现本文算法，并借助可视化工具包(Visualization ToolKit,VTK)^[17] 显示重建结果。分别用α-shape算法和本文算法进行曲面重建，重建结果如图 5所示。其中α-shape算法重建曲面为经过人工多次选择α值获得的相对较好的重建结果。

从图 5(a)可以看出，经过本文算法对原始点云进行简化后，在细节特征丰富的地方点云能得到大量保留，在细节特征较少的地方冗余点云数据得到有效的简化。本文算法根据局部特征尺寸对点云数据进行简化是合理正确的，同时，点云简化步骤使得本文算法能够处理更大数量的点云数据。

如图 5(b)和(c)所示，α-shape算法重建结果中，bunny的后背存在孔洞而在脖子、尾巴和耳朵等地方又有很多冗余的三角面片。blade的右边缘和中间三个横柱处存在多余三角面片，在下面支座的平坦面上存在较多孔洞。dragon的背部凸起和牙齿处有很多冗余三角面片，在身体上存在少量孔洞。而本文算法在这些地方的重建结果能够更加准确地反映物体表面的真实情况。对于非均匀的点云数据，α-shape算法进行曲面重建时，无法使用一个固定的α值来解决重建结果中冗余三角面片和孔洞的矛盾，而本文算法只需要输入原始点云数据，就可以完全自适应地重建出物体表面，重建的表面与物体表面基本一致，且表面几乎没有孔洞。

从表 1可以看出，本文算法能有效减少点云数量，点云简化率能达到70%左右，极大地提高了后续处理步骤的效率。本文算法重建结果中的三角面片数量比α-shape算法相对较少，因此冗余的三角面片也更少。本文算法重建耗时主要由Delaunay三角剖分耗时、点云简化耗时和曲面重建耗时三部分组成，且与点云数量大致呈线性关系。点云简化需要计算所有原始数据点的局部特征尺寸，而且涉及到大量邻近点和球形邻域搜索，故耗费的时间是最多的，曲面重建时无需重新计算点云的局部特征尺寸且只有简化后的点云参与曲面重建，故耗时较少。而α-shape算法的重建耗时是选择n个固定的α值进行n次重建的时间总和，因此α-shape算法重建耗时具有不确定性。

当点云按曲率自适应简化后，文献^[10]算法的重建结果会出现大量而明显的孔洞。为说明本文算法适用性强，采用在geomagic2013中依据曲率简化，简化率为70%的简化点云数据horse为实验数据。因为此数据已经简化，所以使用本文算法时不再对点云进行简化，本文算法的重建结果如图 6所示。

如图 6所示，horse的重建结果与其原始表面基本一致，horse的耳朵和马蹄等区域的重建结果能够反映原始数据的细节特征信息，在腹部和脖子等区域的重建结果含有个别的孔洞，这些孔洞在后续的处理中很容易修补。故本文算法仍然适用于根据曲率简化的非均匀点云的曲面重建，适用范围比文献^[10]算法更广。

4 结语

本文对非均匀散乱点云数据表面重建进行研究，根据采样点的局部特征尺寸改进了α-shape曲面重建算法。算法先根据局部特征尺寸简化点云数据，再自适应地重建物体表面。实验结果表明本文算法能够有效合理准确地减少点云数据，解决了α-shape算法不适用于非均匀点云数据曲面重建的问题，且本文算法适用范围比文献^[10]算法更广。但是本文算法重建出的表面是非流形表面，为了获得更好的重建结果，通过法向过滤三角面片、2-流形面提取和曲面法向一致化等后续处理步骤需要进一步研究。