層次聚類法在逆向工程中的應用研究

白代萍 ， 李艷芹 ， 丁 靜 ， 楊 勇 ， 樊 躍

（1. 鄭州輕工業學院機電工程學院，河南 鄭州 450002；2. 鄭州鐵路局，河南 鄭州 450053）

隨著數控、計算機、三維光學檢測技術的飛速發展，逆向工程中的產品數字化測量手段和方法變得越來越豐富，測量過程也變得越來越簡單。龐大的測量數據（大規模的離散測量點，點云）為描述原型產品的基本形狀特征和結構細節提供了充足的信息[1]。在測量時由于受工作范圍和被測件本身復雜度的影響，將多視圖下的數據轉化為同一坐標系后往往存在大量冗余點，為存儲和后續處理帶來很大不便。因此，在滿足精度的前提下對掃描數據進行采樣是逆向處理和質量在線檢測的重要工作。

國內外學者對逆向工程中的點云采樣技術進行了大量研究。LEE等[2]提出了三維柵格簡化算法，利用八叉樹剖分生成三維柵格，遍歷柵格若其中點云的法向量偏差大于指定閾值則細分單元格，對每個柵格選擇代表點形成簡化點云。PAULY等[3]提出了層次聚類法，自上而下對點云進行二叉分割成多個子集，若子集滿足曲面變分和點的數目大于給定閾值則進一步分割。Alexa等[4]根據點對其局部構建的移動最小二乘曲面影響程度大小，去除影響較小的點來減少冗余。算法能較好地保留重建特征，但涉及非線性求解問題，容易出錯。Song等[5]提出全局聚類采樣，先進行隨機聚類，再對每個類使用搜索算法進行優化，選取合理的特征點使得目標函數值最小以實現點云的簡化。王仁芳等[6]提出了一種基于相似性的自適應點模型簡化算法，將點云劃分成特征部分和非特征部分并分別處理得到簡化點云。杜小燕等[7]提出了一種較新的基于幾何圖像的點云數據簡化算法，將點集合的笛卡爾坐標首先轉換為球面極坐標，再將球面極坐標重采樣到灰度圖像中進行幾何圖像的簡化，圖像灰度值的變化反映了采樣表面曲率變化信息。曹巨明等[8]提出了針對線掃描點云的移動最小二乘增量式多視點云數據融合算法，通過尋找當前點云與已增量式融合的點云數據的重疊部分，在重疊部分數據集上構建移動最小二乘曲面，將重疊部分的每一個在移動最小曲面上的對應點合并到當前以增量式融合的點云集中，從而實現了點云數據的融合。梁新合等[9]采用與法向局部方差有關的自適應角度閾值的截斷函數限制鄰域點的選擇，獲得與表面特征有關的自適應最優鄰域，采用改進的三邊濾波方法實現法向矢量濾波和位置濾波。實驗驗證了該方法的可行性，與其他濾波方法相比，該算法能更有效地保持細節特征，同時獲得光順的離散表面。

本文在深入研究層次聚類法的基礎上，對其進行了分類定義，研究了曲面變分、BSP構建等關鍵技術，構建了算法實現的流程圖。應用該算法結合Geomagic Qualify軟件對某汽輪機廠的大型水輪機葉片進行采樣和逆向建模研究。

1 層次聚類法

1.1 層次聚類法的分類定義

從系統樹圖形成的方式來看，層次聚類算法包括凝聚式算法和分裂式算法。

凝聚式算法以“自底向上”的方式進行。首先將每個樣本作為一個聚類，然后合并相似性最大的聚類為一個新的聚類，直到所有的聚類都被融合成一個新的聚類。它以若干個聚類開始，以1個聚類結束。

分裂式算法以“自頂向下”的方式進行，首先，將整個樣本集合作為一個大的聚類；其次，在算法進行的過程中考察所有可能的分裂方法，把整個聚類分成若干個小的聚類。第1步分成2類，第2步分成3類，依此規律一直進行下去直到最后一步分成n類。在每一步中選擇一個使得相異程度最小的分裂。該方法得到一個相反結構的系統樹圖，以1個聚類開始，以n個聚類結束，與凝聚式算法相比，由于凝聚式算法在計算上簡單、快捷，最終結果相近，所以絕大多數層次聚類方法都是凝聚式的，區別在于聚類的相似性度量的定義有所不同。

層次聚類算法[10]（圖1）的目標是針對具有n個樣本的集合 X ∈ Rn×d，通過相似性函數計算樣本間的相似性并構成相似性矩陣 R = (rij)n×n；再根據樣本間的相似性矩陣把樣本集組成一個分層結構，產生一個從1到n的聚類序列。該序列具有二叉樹的形式，即每個樹的結點有兩個分支，從而使得聚類結果構成樣本集的系統樹圖。

圖1 全局聚類采樣

聚類法廣泛地應用于計算機圖形學中三維復雜實體處理[3]。分裂式算法將點云P分割成大量的子集{Ci}，即 P ={Ci}，每一個子集Ci用一個點pi來取代，形成一個簡化點集 P ′= { pi} 。以一平面點集示例，如圖2所示。線寬對應著分割層次，對點云進行3次劃分，形成 (23-1)=7個聚類，對類 Ci(i ∈ [1 ,7])用點pi取代，生成一個簡化點集由7個點組成取代原來的19個點。

圖2 層次聚類法二維示意圖

1.2 關鍵技術研究

1.2.1 曲面變分

曲面變分是點云局部面的幾何屬性的度量，該參數在層次聚類算法中作為二叉劃分的依據，其計算使用協方差分析[11]求解，它依賴于點的鄰域。考慮點云為無規則的散亂點3，采樣點p∈P的k鄰域記為Np，如圖3所示。

圖3 點云局部面協方差分析二維示意圖

給定采樣點p對應3×3協方差矩陣，

T：表示矩陣的轉置。

考慮對稱半正定協方差矩陣C的特征值問題：

式中，λl：矩陣C的特征值。

特征值λl是實值，評價了pi沿著特征向量方向的偏移。總的偏移也就是鄰域點pi到質心的歐式距離的平方和，如式(3)：

在求實對稱矩陣的特征根與特征向量時，采用雅可比法。雅克比算法的理論基礎是若矩陣A與 B相似，則它們具有相同的特征值；如果 A是實對稱矩陣，則一定可用正交相似變化使其對角化。在程序中，用旋轉矩陣T不斷對矩陣A作正交相似變換把 A化為對角矩陣，從而求出 A的特征值與特征向量。

假定λ0≤λ1≤λ2，定義：

式(4)中，定義平面T(x)過點p且滿足鄰域點到面的距離平方和最小，因此v0逼近面在點p處的法矢np，也就是說v1和v2橫跨點 p處平面T(x)。λ0定量描述了點沿法矢的偏移，λ1和λ2表示了沿切平面主軸方向偏移。

定義點p處的曲面變分：

1.2.2 BSP鄰域構建

給定點云P，對于其中某點pi鄰域可以分為多種。本文著重介紹BSP鄰域，它的計算使用到二叉樹空間分割。

該算法認為：

任何平面都可以被分割成兩個半空間。位于平面一側的點定義了一個半空間，位于另一側的點定義了另一個半空間。在半空間中又存在一個平面，將此半空間進一步分割成兩個更小的子空間。依次循環，得到越來越小的子空間，這些空間構造成一個二叉樹[12]。

在二叉樹中，初始點云對應樹的根結點，位于同一子空間中的點云對應某一葉子結點。結點在分割時必須滿足以下兩個條件：

（1）結點包含點的數目大于用戶指定類所能包含最大點數nmax；

（2）結點的曲面變分σ (P)大于給定曲面變分閾值σmax。

分割平面由點云質心和 P對應協方差矩陣的特征向量v2來確定，特征向量v2對應協方差矩陣的最大特征根，使得點云總是沿著變化最大的方向分割。不滿足分割判據時，點云將形成一個葉子結點對應一個類Ci[13]。

2 算法描述實現

在1.1節中已經提到，層次聚類是一種自上向下的劃分點云方法，實際上是空間二叉分割，是一遞歸過程。空間劃分結束時，每個葉子結點對應一個類，類的總數等于采樣后的點數。每個類用均值點取代，得到簡化點云。算法流程說明如圖4所示。

圖4 層次聚類法算法流程

（1）讀入數據，初始化；

（2）建立結點的協方差矩陣，求解結點的曲面變分；

（3）判斷結點是否滿足分割條件，若不滿足則生成類即葉子結點，并跳至步驟(6)；

（4）計算結點分割方向，判斷結點中各點屬于左或者右子樹；

（5）左、右子樹形成新的結點；

（6）判斷是否遍歷所有結點，若有未遍歷結點轉至步驟(2)；

（7）遍歷各類即葉子結點，一般用類質心取代該類；

（8）采樣結束。

對圖5(a)所示的兔子點云（34384個點）進行聚類劃分，一個顏色區域代表最終顯示結果為一個點，點位于區域的中心，簡化為圖5(b)（1000個點），從該圖看出其仍然顯著保持了原幾何物體的特征信息。

圖5 聚類云圖和點云顯示[11]

3 實 驗

該項目是對某汽輪機廠生產的水輪機葉片進行數字化質量檢測，傳統的葉片檢測主要采用三坐標測量和標準樣板法。三坐標測量精度較高，但每次只能完成一個測點，測量葉片耗時較長，制約了葉片生產效率。標準樣板法需要針對不同級的葉片對應不同檢驗模板數目，這種方法費時費事，而且引入較大測量誤差。光學設備的發展為葉片的數字化檢測提供了基礎，主要原理是：通過光學設備掃描到葉片的點云數據，對這些數據進行預處理，處理完畢的數據與工件模型進行比對。流程如圖6所示。其中圖6 (a)是采集到的原始點云，圖6 (b)是使用層次聚類法，采樣后的點云。圖6 (c)是在Geomagic Qualify比對軟件中，對偏差進行評價用色譜圖顯示。

圖6 水輪機葉片數字化檢測

葉片模型采集到的原始點云有180萬點，使用層次聚類法 1/70采樣后得到的點云為 2.6萬點，和傳統三坐標機只能打出最多幾百個點相比，能滿足測點數量的要求。計算簡化點云和原始點云的平均偏差0.746 mm，數據和模型匹配時采用ICP精確配準，并正確建立點云與CAD模型關系，相比傳統人工方法，具有較高精度，能滿足檢測要求。

4 結 論

（1）層次聚類采樣算法對點云自上而下進行二叉劃分，劃分結束用均值點來取代各類，得到簡化點云。二叉劃分依據是曲面變分和類中點的數目，曲面變分的計算是根據點的鄰域，構建3×3協方差矩陣，使用雅克比迭代法求解矩陣特征值和特征根。

（2）層次聚類采樣算法的優點在于保留點較均勻，計算速度較快，但構建協方差矩陣時至少要有4個點，也就是二叉劃分形成的類要包括3個以上點，不能實現任意比率采樣，適用于大比率保留點較少的場合，因此需要進一步研究其他采樣算法。

（3）采用層次聚類采樣算法對東風汽輪機廠生產的水輪機葉片進行數字化質量檢測，點云關系和CAD數模對比證明了該采樣算法的正確性，可以應用在逆向工程上。

[1] 柯映林. 反求工程CAD建模理論、方法和系統[M].北京: 機械工業出版社,2005: 16-17.

[2] Lee K H,Woo H,Suk T. Point data reduction using 3D grids [J]. The International Journal of Advanced Manufacturing Technology,2001,18 (3): 201-210.

[3] Pauly M,Gross M,Kobbelt L P. Efficient simplification of point-sampled surfaces[C]//IEEE Visualization,2002,2002: 163-170.

[4] Alexa M,Behr J,Cohen-Or D,Fleishman S,Levin D,Silva C T. Point set surfaces[C]//Visualization,2001.VIS′01. Proceedings,IEEE of Computer Society,2001: 21-29,537.

[5] Song Hao,Feng H Y. A global clustering approach to point cloud simplification with a specified data reduction ratio [J]. Computer-Aided Design,2008,40(3): 281-292.

[6] 王仁芳,張三元,葉修梓. 基于相似性的點模型簡化算法[J]. 浙江大學學報(工學版),2009,43(3):448-454.

[7] 杜小燕,郝傳剛,王玉梅. 基于幾何圖像的點云數據簡化算法[J]. 蘇州大學學報(工科版),2009,29(2): 13-16.

[8] 曹巨明,吾守爾·斯拉木,梁 晉,梁新合,張德海.移動最小二乘增量式多視點云數據融合算法[J]. 西安交通大學學報,2009,43(9): 46-50.

[9] 梁新合,梁 晉,郭 成,曹巨明. 基于自適應最優領域的散亂點云降噪技術研究[J]. 中國機械工程,2010,21(6): 639-643.

[10] Sambasivam S,Theodosopoulos N. Advanced data clustering methods of mining web documents [J].Issues in Informing Science and Information Technology,2006,8(3): 563-579.

[11] Pfister H,Zwicker M,Van Baar J,Gross M. Surfels:Surface elements as rendering primitives [C]//ACM Press/Addison-Wesley Publishing Company,2000:335-342.

[12] Bentley J L. Multidimensional binary search trees used for associative searching [J]. Communications of the ACM,1975,18(9): 509-517.

[13] Shamir A,Shapira L,Cohen-Or D. Mesh analysis using geodesic mean-shift [J]. The Visual Computer,2006,22(2): 99-108.