局部相關(guān)點(diǎn)平衡的三維激光點(diǎn)云骨架提取

周明全，李純輝，王麗清，張雨禾，耿國(guó)華

（西北大學(xué) 信息科學(xué)與技術(shù)學(xué)院，陜西 西安 710127）

1 引言

通過(guò)激光掃描設(shè)備掃描獲取物體的點(diǎn)云模型，具有簡(jiǎn)單便捷、準(zhǔn)確率高等優(yōu)點(diǎn)［1］，成為目前最常用的數(shù)字化手段之一［2］。由一組骨架點(diǎn)和骨架點(diǎn)之間的拓?fù)潢P(guān)系［2-4］構(gòu)成的點(diǎn)云模型骨架線能夠有效描述點(diǎn)云模型的整體形狀和拓?fù)浣Y(jié)構(gòu)［4-5］。因此，骨架線的計(jì)算、提取和匹配等成為計(jì)算機(jī)圖形學(xué)［6-7］、計(jì)算機(jī)視覺(jué)［8-9］和機(jī)器智能［10］等領(lǐng)域的研究熱點(diǎn)，并被廣泛應(yīng)用于形狀分析［11］與重建［12］、路徑規(guī)劃［13-14］、動(dòng)作和手勢(shì)識(shí)別［15-17］、模型檢索［18-20］以及植物分析［21-22］等任務(wù)中。基于骨架分析的相關(guān)研究都需要以骨架提取為基礎(chǔ)，因此骨架線的提取至關(guān)重要。常見(jiàn)的骨架類(lèi)型按照空間維度可以分為2D骨架和3D骨架兩類(lèi)［12］，其中2D骨架是由一組位于模型各分支的中間曲 線 構(gòu) 成 的［4，6，23］，而3D骨架［24］則是由 中 心曲面和中心曲線共同構(gòu)成的。相較于3D骨架，2D骨架線具有信息壓縮程度高、拓?fù)潢P(guān)系明確等優(yōu)勢(shì)，因此其應(yīng)用也更為廣泛［9］。

目前已有諸多研究人員開(kāi)展了點(diǎn)云模型2D骨架線提取相關(guān)工作，具體可以分為：（1）基于質(zhì)心Voronoi細(xì) 分（Centroidal Voronoi Tessellation，CVT）的方法［25-26］。該方法是骨架提取中應(yīng)用最多的策略。例如：Lu等［26］通過(guò)變分特征將CVT概念推廣到圖中，將CVT應(yīng)用在骨架提取的任務(wù)中，在提取骨架的同時(shí)還實(shí)現(xiàn)了對(duì)模型的形狀分割，但該方法時(shí)間復(fù)雜度較高，針對(duì)大規(guī)模點(diǎn)云模型處理時(shí)間較長(zhǎng)。（2）基于平均曲率流（Mean Curvature Flow，MCF）的方法［27-29］，該方法能提取到高質(zhì)量的骨架。例如：Chuang等［29］發(fā)表了一種基于有限元跟蹤的快速M(fèi)CF的方法，隨 后Tagliasacchi等［27］將MCF應(yīng) 用 在 骨 架 提 取的任務(wù)中，通過(guò)分析流動(dòng)的微分特征，發(fā)現(xiàn)MCF局部增加了形狀各向異性，利用其面積最小化的特性將曲率流推向極值，從而折疊輸入網(wǎng)格幾何體并獲得骨架結(jié)構(gòu)。（3）基于拉普拉斯收縮的算法。Cao等［4］提出了一種基于拉普拉斯收縮的曲線骨架提取算法（Laplacian-Based Contraction，LBC）。但是這些方法只能應(yīng)用在經(jīng)過(guò)降噪、去除異常值等操作處理之后的模型上，難以應(yīng)用于使用激光掃描設(shè)備所獲取的點(diǎn)云模型。（4）基于深度學(xué)習(xí)的方法。Atienza等［30］基于Pyramid UNets提出了PSPU-SkelNet，適用于從二維形狀點(diǎn)云模型中提取曲線骨架任務(wù)。Lin等［24］基于中軸變換（Medial Axis Transform，MAT）提出了Point2Skeleton（P2S），用于從三維點(diǎn)云模型中提取二維骨架，從而指導(dǎo)表面重構(gòu)和分割。Dou等［31］基于覆蓋軸（Coverage Axis）提出一種深度學(xué)習(xí)的提取骨架的算法，其結(jié)果接近基于MAT的骨架提取算法。

適用于激光掃描點(diǎn)云模型的曲線骨架提取方法，大多都基于幾何學(xué)相關(guān)原理［2-7］。Tagliasacchi等［2］基于旋轉(zhuǎn)對(duì)稱(chēng)軸（Rotational Symmetry Axis，ROSA）的概念，構(gòu)造了一種可以針對(duì)缺失大量數(shù)據(jù)的不完整點(diǎn)云模型的曲線骨架提取方法。Huang等［6］將L1中軸（L1-Medial）應(yīng)用于曲線骨架提取任務(wù)中，將L1中軸局部化自適應(yīng)提取不同位置的中值點(diǎn)作為骨架點(diǎn)。Qin等［7］通過(guò)構(gòu)建基于質(zhì)量驅(qū)動(dòng)拓?fù)涓兄淖顑?yōu)傳輸方案提取粗糙曲線骨架（MdCS），然后使用能量正則化完成幾何優(yōu)化，最終獲得平滑曲線骨架。ROSA雖然能夠應(yīng)對(duì)大面積缺失數(shù)據(jù)的原始掃描點(diǎn)云模型，但是旋轉(zhuǎn)對(duì)稱(chēng)軸的計(jì)算需要無(wú)噪聲擾動(dòng)的點(diǎn)云，因此針對(duì)噪聲點(diǎn)云效果較差。L1中軸骨架和MdCS雖然能應(yīng)對(duì)存在噪聲、異常值和缺失數(shù)據(jù)的原始掃描點(diǎn)云模型，但是L1中軸骨架不能處理多分支結(jié)構(gòu)的模型，MdCS拓?fù)浣Y(jié)構(gòu)準(zhǔn)確度不高，且計(jì)算所需時(shí)間資源消耗巨大。

針對(duì)上述工作中存在的問(wèn)題，受文獻(xiàn)［32-34］的啟發(fā)，本文提出一種快速提取激光掃描點(diǎn)云骨架線的方法。該方法將初始骨架點(diǎn)定義為三維模型中具有對(duì)稱(chēng)法向的最近相關(guān)點(diǎn)的中點(diǎn)，通過(guò)中點(diǎn)狀態(tài)的平衡化迭代得到最終的骨架點(diǎn)；然后使用主成分分析法（Principal Component Analysis，PCA）分析骨架點(diǎn)的方向一致性，從而得到骨架曲線段；最后將使用廣度優(yōu)先搜索（Breadth First Search，BFS）合并后的曲線段利用拉普拉斯平滑進(jìn)行優(yōu)化并得到最終曲線骨架。本文實(shí)驗(yàn)將此算法應(yīng)用在通過(guò)激光掃描設(shè)備得到的點(diǎn)云模型上，驗(yàn)證了算法的有效性、魯棒性和高效性。將提取的激光掃描點(diǎn)云的骨架線應(yīng)用在點(diǎn)云模型的動(dòng)作變換任務(wù)中，展示了骨架線的實(shí)際作用。

2 基本原理

本文的骨架提取方法主要可以分為兩步：提取骨架點(diǎn)和生成骨架線。首先，通過(guò)選擇具有對(duì)稱(chēng)法向和最近距離的對(duì)稱(chēng)點(diǎn)的中點(diǎn)，獲取初始骨架點(diǎn)，如圖1（a）所示；應(yīng)用射線法去除模型外的初始骨架點(diǎn)，如圖1（b）所示；通過(guò)迭代來(lái)平衡骨架點(diǎn)并篩選權(quán)重大的骨架點(diǎn)組成最終骨架點(diǎn)，如圖1（c）所示。然后，基于PCA判斷骨架點(diǎn)的方向一致性并重建骨架曲線，采用廣度優(yōu)先遍歷算法合并骨架曲線，效果如圖1（d）所示；最后將平滑后的骨架曲線相連得到最終的曲線骨架，如圖1（e）所示；基于骨架的模型動(dòng)作變換實(shí)際效果如圖1（f）所示。

圖1 本文曲線骨架提取算法的過(guò)程展示Fig.1 Pipeline of proposed method

2.1 提取骨架點(diǎn)

骨架點(diǎn)可以看作模型的部分采樣點(diǎn)投影到模型內(nèi)部得到的局部中心點(diǎn)［6］，也就是說(shuō)，初始的骨架點(diǎn)是模型的點(diǎn)集通過(guò)采樣處理得到的子集；另外，直接對(duì)采樣點(diǎn)進(jìn)行計(jì)算，還可以有效提高骨架線提取的效率，避免遍歷點(diǎn)云上的所有點(diǎn)數(shù)據(jù)。除了基于效率的考慮，對(duì)模型下采樣還能減少模型本身的一些不利因素對(duì)最終骨架的影響。例如，當(dāng)模型不均勻時(shí)，通過(guò)均勻采樣便可以避免這一缺陷；當(dāng)遇到噪聲時(shí)，采樣也會(huì)減少異常值噪聲對(duì)骨架提取任務(wù)的干擾。因此，本文算法首先進(jìn)行隨機(jī)采樣，使得采樣點(diǎn)能均勻的覆蓋模型，從而避免骨架線變形。具體而言，對(duì)于一個(gè)離 散 點(diǎn)云模型P={pi}i∈I?R3，其中I是 點(diǎn) 云P的點(diǎn)數(shù)，將采樣點(diǎn)組成的集合記作Ps，從P中使用隨機(jī)算法獲得點(diǎn)的索引，本文實(shí)驗(yàn)中采樣點(diǎn)的規(guī)模為0.1×I。

2.1.1 法向估計(jì)和重定向

本文中初始骨架點(diǎn)需要基于法向求解，所以首先需要估計(jì)采樣點(diǎn)的法向［35-36］。本文使用PCA對(duì)Ps中的點(diǎn)進(jìn)行法向估計(jì)，即：ps及其在P上的鄰域點(diǎn)共同構(gòu)成點(diǎn)集Pt，對(duì)Pt進(jìn)行特征分解得到特征值和特征向量，將最小特征值對(duì)應(yīng)的特征向量pn t作為ps的法向。具體如下，利用式（1）計(jì)算鄰域點(diǎn)集Pt的中心點(diǎn)pm：

此時(shí)ps和法向pnt應(yīng)該滿(mǎn)足pnt使 得Pt中 所 有的鄰域點(diǎn)在方向pnt上的投影分布最密集：

式中：k=20，將pt-pm記作則有

由于上述過(guò)程計(jì)算的法向具有二義性［37］，如圖2（a）所示，而這種二義性會(huì)影響骨架點(diǎn)的求解過(guò)程，所以還需要對(duì)法向的二義性進(jìn)行消除，即法向的重定向。具體而言，對(duì)于兩個(gè)距離很近的采樣點(diǎn)pi，pj，如果他們的法向pni，pnj同向，則滿(mǎn)足≈1。本文使用深度優(yōu)先算法（Depth First Search，DFS）來(lái)快速完成重定向。具體而言，選擇距離Pt的中心點(diǎn)PsO最遠(yuǎn)的采樣點(diǎn)pi，如果‖(pi+pni)-PsO‖＞‖pi-PsO‖，則pni=-pni，反 之則不改變pni。隨后使用K最鄰近算法（K-Nearest Neighbor，KNN）可以直接得到pi的距離最近的一個(gè)Ps中的點(diǎn)pj，如果pni?pnj/(|pni|×|pnj|)≤-0.5，則pnj=-pnj，并將pj標(biāo)記為已經(jīng)訪問(wèn)的狀態(tài)，此時(shí)令pi=pj。通過(guò)KNN算法尋找一個(gè)沒(méi)有被標(biāo)記為已訪問(wèn)狀態(tài)的采樣點(diǎn)作為pj，重復(fù)上述過(guò)程，直到Pt中所有的點(diǎn)都被標(biāo)記為訪問(wèn)的狀態(tài)。此時(shí)Pt中所有的點(diǎn)都將指向模型的同側(cè)，消除了法向的二義性，如圖2（b）所示。用該重定向方法重定向之后的法向，都指向模型的內(nèi)部。

圖2 法向二義性示意圖Fig.2 Schematic diagram of normal vector ambiguity

2.1.2 骨架點(diǎn)的求解

要進(jìn)行骨架點(diǎn)識(shí)別，本節(jié)首先做以下定義：

定義一：對(duì)于采樣點(diǎn)集合Ps中的任意一點(diǎn)pi，存在pj∈P，j∈S{i}是和pi具有相反法向的點(diǎn)中距離pi最近的點(diǎn)，將pj和pi定義為一組相關(guān)點(diǎn)。

定義二：一組相關(guān)點(diǎn)的中點(diǎn)是初始骨架點(diǎn)。

對(duì)于Ps中的任意一點(diǎn)pi，從Ps中識(shí)別出一組法向與pni相反的點(diǎn)，在這里法向pni和pnj使用一個(gè)帶閾值的點(diǎn)積來(lái)判斷。在得到的這一組法向中，進(jìn)一步找到特定點(diǎn)pj，使得pj到pi的歐氏距離最小，隨后計(jì)算出pj和pi的中點(diǎn)作為初始骨架點(diǎn)。當(dāng)模型的分支之間的距離較小時(shí)，識(shí)別到的相關(guān)點(diǎn)可能并不在一個(gè)分支上，導(dǎo)致提取到無(wú)效的冗余骨架點(diǎn)，如圖1（a）中手臂和身體之間產(chǎn)生的無(wú)效骨架點(diǎn)，此時(shí)需要通過(guò)算法將不在模型中的相關(guān)點(diǎn)之間的中點(diǎn)，也就是無(wú)效的骨架點(diǎn)去除掉。具體而言，本文使用射線法來(lái)去除無(wú)效骨架點(diǎn)，從而得到有效骨架點(diǎn)。射線法具體是指從任意一點(diǎn)出發(fā)，從該點(diǎn)射出的射線，如果與模型邊界的交點(diǎn)數(shù)量為奇數(shù)，則該點(diǎn)在模型內(nèi)部；反之，如果交點(diǎn)為偶數(shù)，則該點(diǎn)在模型的外部，如圖3所示（彩圖見(jiàn)期刊電子版）。

圖3 （a）中所需判斷的點(diǎn)（被標(biāo)記為紅色）在模型外部，黃色射線與模型的交點(diǎn)（被標(biāo)記為綠色）有兩處，是偶數(shù)。（b）中所需判斷的點(diǎn)（被標(biāo)記為紅色）處于模型內(nèi)部，黃色射線與模型的交點(diǎn)（被標(biāo)記為綠色）有一處，為奇數(shù)Fig.3 Figure（a）shows that when the point（marked red）to be judged is outside the model，the intersection of the yellow ray and the model（marked green）has two even numbers.The point（marked red）to be judged in the figure（b）is inside the model，and there is an odd number at the intersection of the yellow ray and the model（marked green）

在將無(wú)效的骨架點(diǎn)去除之后，效果如圖1（b）所示。當(dāng)模型足夠簡(jiǎn)單時(shí)，例如圓柱模型，這些點(diǎn)可以直接形成模型P的骨架，但是當(dāng)模型不是簡(jiǎn)單的基本幾何形狀時(shí)，還需要使初始骨架點(diǎn)進(jìn)一步收斂到最終骨架點(diǎn)的位置。由于此時(shí)的有效骨架點(diǎn)qi∈Q可以看作是模型上的采樣點(diǎn)pi經(jīng)過(guò)位置變換得到的，受文獻(xiàn)［38］的啟發(fā)，在垂直于骨架點(diǎn)法向qni的平面上，在每個(gè)有效的骨架點(diǎn)qi處均勻地發(fā)射出m條射線，其中第j條和第j+m/2條射線的方向是相反的。從骨架點(diǎn)qi出發(fā)，沿著每條射線方向迭代地向前搜索，直到靠近模型上的任意一點(diǎn)，將qi到該點(diǎn)的距離記作dj，將qi到模型上該點(diǎn)的方向標(biāo)記為向量vqji。本文實(shí)驗(yàn)中，m的數(shù)量為16條，被分為8對(duì)，每對(duì)由兩個(gè)方向相反的向量組成。對(duì)于點(diǎn)集Q中的任意一點(diǎn)qi，定義qi的能量Eqi為：

由上述定義可得到骨架點(diǎn)Q的平衡狀態(tài)如式（4）所示：

其中：NQ是Q中點(diǎn)的數(shù)量，r是模型的最小邊界球的半徑。根據(jù)上述定義，可知qi和qi+1之間的關(guān)系如式（5）所示：

其中，L=0.01是步長(zhǎng)參數(shù)。當(dāng)式（3）中能量梯度為0時(shí)，意味著任意一個(gè)qi的E=0，也就是qi的位置處于模型橫截面的中心位置，如圖4（c）所示。根據(jù)上述過(guò)程，便可以完成骨架點(diǎn)的提取。在獲取初始骨架點(diǎn)下個(gè)位置的同時(shí)也可以得到下次移動(dòng)的方向，可以加速初始骨架點(diǎn)平衡化的過(guò)程。但是此時(shí)的骨架點(diǎn)數(shù)量還是遠(yuǎn)遠(yuǎn)大于最終用來(lái)生成骨架點(diǎn)的數(shù)量，所以還需要進(jìn)一步從骨架點(diǎn)中篩選合適的骨架點(diǎn)。這里使用二次鄰域中權(quán)值比較的方式來(lái)完成骨架點(diǎn)的選擇。具體而言，對(duì)于一個(gè)骨架點(diǎn)qi，其權(quán)值被定義為骨架點(diǎn)qi的rknn鄰域內(nèi)點(diǎn)的數(shù)量N。對(duì)于Q中任意一個(gè)骨架點(diǎn)qi，將qi的rknn鄰域中的點(diǎn)記作Q r，則對(duì)于Q r中所有的骨架點(diǎn)而言，應(yīng)該保留的qr應(yīng)該滿(mǎn)足式（6）：

圖4 （a）展示的是得到垂直于qni的平面的方法，pi+1是pi距離最近的采樣點(diǎn)，（b）中骨架點(diǎn)qi處于不平衡的位置，（c）中骨架點(diǎn)qi處于平衡位置Fig.4（a）shows the method of obtaining the plane perpendicular to qni，and pi+1 is the nearest sampling point from pi，in（b），the skeleton point qi is in an unbalanced position，and in（c），the skeleton point qi is in an balanced position

2.2 骨架生成

上述迭代過(guò)程能夠得到一組骨架點(diǎn)Q={qi}，i∈NQ，它是輸入點(diǎn)云模型P的各個(gè)分支的中間地帶，僅表示了模型P骨架的位置信息，而完整的骨架線還包含了骨架點(diǎn)之間及骨架分支之間的拓?fù)溥B接信息。所以本節(jié)計(jì)算骨架點(diǎn)Q的拓?fù)溥B接信息。

首先，使用PCA來(lái)檢查qi和它的鄰域點(diǎn)集是否滿(mǎn)足方向一致性的條件，隨后將滿(mǎn)足要求的鄰域點(diǎn)集生成模型曲線骨架拓?fù)潢P(guān)系的子集，進(jìn)一步通過(guò)廣度優(yōu)先搜索算法將存在連接關(guān)系的子集合并得到模型的骨架曲線分支，再通過(guò)拉普拉斯平滑骨架曲線分支，最后將平滑后的骨架曲線分支組織成完整的曲線骨架。

PCA常被用來(lái)檢測(cè)骨架點(diǎn)是否達(dá)到收斂要求［6，23］，即：對(duì)于骨架點(diǎn)qi，計(jì)算qi及其k鄰域內(nèi)點(diǎn)集的協(xié)方差矩陣的特征值和特征向量，如式（7）所示：

其中：所有的特征值滿(mǎn)足λ0i≤λ1i≤λ2i都是實(shí)值，對(duì)應(yīng)的特征向量{v0i，v1i，v2i}形成一個(gè)正交框架，即qi和其鄰域點(diǎn)的主成分。定義qi和其鄰域點(diǎn)的方向一致性為σi，則有如式（8）的關(guān)系：

σi的實(shí)際含義是當(dāng)σi越接近1，即λ0i，λ1i和λ2i相比越小，更多的點(diǎn)分布在沿著分支對(duì)齊的位置，也就是方向一致性越高。實(shí)踐中，設(shè)置k=3，對(duì)于任意qi而言，如果σi＞0.9，便可以認(rèn)為qi和它的兩個(gè)鄰域點(diǎn)已經(jīng)滿(mǎn)足骨架點(diǎn)的收斂要 求，可 以 認(rèn) 為 對(duì) 于 點(diǎn)qi和qi+1，qi+1和qi+2而言，存在拓?fù)潢P(guān)系。隨后遍歷所有骨架點(diǎn)，便可以得到達(dá)到要求的所有的兩個(gè)點(diǎn)之間的拓?fù)潢P(guān)系，記 作R={(rmn)}m∈Q，n∈Q，m≠n，rmn=(qm，qn)，如圖5（a）所示。

圖5 R和RL的區(qū)別，（a）是R，（b）是根據(jù)R通過(guò)廣度優(yōu)先搜索得到RLFig.5 Difference between R and RL.（a）is R，and（b）is RLobtained through breadth first search according to R

得到R之后，使用廣度優(yōu)先搜索，增長(zhǎng)(qm，qn)的關(guān)系，得到RL，如圖5（b）所示。增長(zhǎng)后的效果如圖1（c）所示。隨后對(duì)于RL中任意rnm=(qm，qm+1，…，qn-1，qn)兩端的端點(diǎn)qm和qn，對(duì)qm尋找最近骨架點(diǎn)qm+1，構(gòu)成方向向量vm+1m=qm+1-qm，同理得qn的方向向量vn n-1=qn

qn-1，使它們?cè)诖朔较蛏系鷮ふ液线m的骨架點(diǎn)相連即可得到完整的骨架。由于在此過(guò)程中，在尺度變化劇烈時(shí)，可能會(huì)導(dǎo)致骨架曲線的變化不夠平滑。為了減少這種劇烈變化帶來(lái)的質(zhì)量上的影響，采用拉普拉斯收縮以獲得平滑的骨架曲線。對(duì)于拉普拉斯平滑策略點(diǎn)qm的定義如式（9）所示：

其中：i=m+1，m+2，…，n-1。

經(jīng)過(guò)平滑后，便可以得到完整的曲線骨架，如圖1（d）所示。

3 實(shí)驗(yàn)以及結(jié)果分析

本文實(shí)驗(yàn)所用硬件和實(shí)驗(yàn)環(huán)境包括：CPU是AMD的 銳 龍7代3700X，內(nèi) 存 為32 GB，GPU為NVIDIA GeForce RTX 2060 SUPER，Visual Studio 2017 Professional，PCL（Point Cloud Library）1.8.0，開(kāi)發(fā)語(yǔ)言為C++。為了測(cè)試和評(píng)估算法的實(shí)際性能和效果，本文實(shí)驗(yàn)采用了公開(kāi)數(shù)據(jù)集Visionair Repository中的多個(gè)三維模型，以及部分利用FARO Focus Premium激光掃描儀激光掃描設(shè)備掃描獲取的點(diǎn)云模型作為測(cè)試數(shù)據(jù)，測(cè)試本文方法及對(duì)比方法的結(jié)果。

3.1 本文方法實(shí)驗(yàn)結(jié)果

圖6所示為本文算法應(yīng)用于鹿、手掌、人、雕塑模型的實(shí)驗(yàn)結(jié)果。圖6（a）中鹿模型的難點(diǎn)在于，鹿的脖子部位存在軸向偏轉(zhuǎn)，骨架點(diǎn)位置易偏移，導(dǎo)致骨架線失真；然而，從圖6中的結(jié)果可以看出本文算法受到的擾動(dòng)較小，原因是雖然存在軸向偏轉(zhuǎn)，但軸向偏轉(zhuǎn)并沒(méi)有導(dǎo)致模型脖子部位的形狀產(chǎn)生較大變化，所以對(duì)骨架生成沒(méi)有產(chǎn)生干擾。圖6（b）中手掌模型的難點(diǎn)在于多個(gè)距離很近的手指會(huì)導(dǎo)致骨架點(diǎn)出現(xiàn)在手指之間的位置，產(chǎn)生錯(cuò)誤的骨架點(diǎn)和骨架線。在本文中，為了避免這種現(xiàn)象的產(chǎn)生，采用了射線法來(lái)去除掉出現(xiàn)在模型之外的錯(cuò)誤位置的骨架點(diǎn)，最終結(jié)果證實(shí)了此算法的有效性。圖6（c）中人模型對(duì)于骨架提取算法的挑戰(zhàn)在于可能會(huì)在手臂和身體之間產(chǎn)生錯(cuò)誤的分支；圖6（d）中雕塑模型的難點(diǎn)在于存在復(fù)雜的拓?fù)浣Y(jié)構(gòu)，然而本文算法都能有效應(yīng)對(duì)這些難點(diǎn)，提取出完整平滑的曲線骨架。

圖6 本文算法提取的模型的曲線骨架（右側(cè)）和原模型（左側(cè)）Fig.6 Curve skeleton of some models（right）and the display of the original model（left）extracted by this algorithm

為了進(jìn)一步驗(yàn)證本文算法的有效性和優(yōu)勢(shì)，本文將其與目前3種前沿方法進(jìn)行比較，即：L1中軸骨架算法、MdCS算法和Point2Skeleton算法。Point2Skeleton所使用的模型已按照其要求通過(guò)均勻采樣，將點(diǎn)云下采樣到2 000點(diǎn)后作為其輸入。圖7（a）為4種算法應(yīng)用于稀疏點(diǎn)云狗模型上的結(jié)果。由于狗模型本身的點(diǎn)數(shù)較少，狗模型的腹部存在一些點(diǎn)的空白，導(dǎo)致使用點(diǎn)的坐標(biāo)作為輸入的方法提取到曲線骨架的難度更高。L1中軸骨架算法由于內(nèi)切球的半徑增長(zhǎng)，導(dǎo)致無(wú)法提取到其他腿部的骨架。MdCS算法在狗模型的身體上提取到了錯(cuò)誤的拓?fù)浣Y(jié)構(gòu)。Point2 Skeleton算法將狗的身體提取為二維骨架，雖然存在拓?fù)溥B接錯(cuò)誤和狗模型的右前腿的骨架缺失，但整體效果好于L1中軸骨架算法。然而，本文算法能夠得到較為滿(mǎn)意的效果，這是因?yàn)殡m然模型上點(diǎn)缺失的部分難以估計(jì)法向，但本文采用了迭代提取骨架點(diǎn)的方法，降低了這種稀疏采樣對(duì)骨架線提取的影響。圖7（b）所示為馬模型的骨架線提取結(jié)果。馬模型在骨架線提取任務(wù)中的難點(diǎn)在于，四肢和身軀之間存在巨大的尺度變化，這種尺度變化會(huì)影響基于幾何性質(zhì)收縮或者內(nèi)切球增長(zhǎng)的算法的有效性。從圖中可以看出，L1中軸骨架算法繼續(xù)受到內(nèi)切球半徑增長(zhǎng)導(dǎo)致分支骨架缺失，MdCS算法即使存在拓?fù)錂z查的過(guò)程，但依然不能應(yīng)對(duì)馬腿和身體之間的巨大尺度變化。Point2Skeleton算法在馬頭部位與真實(shí)情況相差較大，在身體部位提取的二維骨架失真嚴(yán)重。而本文算法能夠得到合理的拓?fù)潢P(guān)系并生成完整的骨架線。圖7（c）所示為人模型，人腿部和腹部的連接點(diǎn)的選擇是一個(gè)難點(diǎn)，導(dǎo)致L1中軸骨架算法選擇的連接點(diǎn)不在模型的內(nèi)部，MdCS算法在腹部依然存在拓?fù)浣Y(jié)構(gòu)錯(cuò)誤，Point2Skeleton算法則未提取到模型的左手。而本文算法在腹部未選擇合并，導(dǎo)致腿部骨架和身體骨架合二為一。

圖7 各種方法提取的曲線骨架效果對(duì)比Fig.7 Effect comparison of curve skeleton extracted by various methods

3.2 魯棒性實(shí)驗(yàn)

為了驗(yàn)證本文算法的魯棒性，本節(jié)針對(duì)添加高斯噪聲的模型開(kāi)展實(shí)驗(yàn)。本節(jié)實(shí)驗(yàn)中，定義高斯噪聲的標(biāo)準(zhǔn)差為2×dmin。其中dmin是模型中點(diǎn)之間最小距離。結(jié)果如圖8所示。MdCS算法在應(yīng)對(duì)有噪聲的點(diǎn)云模型時(shí)效果好于沒(méi)有噪聲的模型，原因在于，添加噪聲后，模型點(diǎn)的坐標(biāo)發(fā)生了變化，導(dǎo)致算法在拓?fù)錂z查時(shí)使用了較大的半徑，減少了拓?fù)溴e(cuò)誤。由于MAT易受噪聲影響，所以Point2Skeleton算法在應(yīng)對(duì)噪聲時(shí)的表現(xiàn)變差，產(chǎn)生錯(cuò)誤拓?fù)溥B接的可能性相較于不受噪聲干擾時(shí)有所提高。受到噪聲干擾時(shí)，Point2Skeleton生成了更少的拓?fù)溥B接，由此雖然拓?fù)溥B接錯(cuò)誤的數(shù)量下降了，但是骨架的質(zhì)量卻沒(méi)有提高。

圖8 輸入帶噪聲的模型時(shí)各種方法提取的曲線骨架效果對(duì)比Fig.8 Effect comparison of curve skeleton extracted by various methods when inputting the model with noise

由于較低程度的高斯噪聲對(duì)PCA中最小特征值的影響有限，所以對(duì)基于PCA的法向計(jì)算的影響也非常有限。并且中點(diǎn)狀態(tài)的平衡化迭代可以在得到居中性更好地骨架點(diǎn)的同時(shí)，降低噪聲的干擾，所以本文算法受到噪聲影響較小。

為了定量地分析幾種算法的實(shí)際效果，統(tǒng)計(jì)了幾種方法在應(yīng)對(duì)模型無(wú)噪聲和有噪聲的情況下的骨架點(diǎn)數(shù)量、拓?fù)溥B接的數(shù)量、交叉點(diǎn)數(shù)量和拓?fù)溥B接錯(cuò)誤的數(shù)量。其中，拓?fù)溥B接錯(cuò)誤的數(shù)量等于骨架缺少的分支數(shù)量和去掉的拓?fù)溥B接數(shù)量的和，去掉的拓?fù)溥B接數(shù)量為去掉拓?fù)淠苁雇負(fù)溥B接正確所需去除的數(shù)量。如表1所示。從表1可以看出，L1中軸骨架線提取算法提取的骨架點(diǎn)數(shù)量較多，拓?fù)溥B接錯(cuò)誤的數(shù)量也較多，骨架分支缺失的情況較多。MdCS算法由于提取的骨架線可能存在環(huán)結(jié)構(gòu)，所以拓?fù)溥B接錯(cuò)誤較多，而且可能存在骨架點(diǎn)和拓?fù)溥B接在模型外的情況。Point2Skeleton具有固定的輸出，即：100個(gè)骨架點(diǎn)和骨架點(diǎn)的拓?fù)潢P(guān)系。由于其輸出的是二維骨架，二維骨架的連接狀態(tài)較多，所以拓?fù)溥B接關(guān)系出現(xiàn)錯(cuò)誤的可能性較大，這也反應(yīng)在了圖8中。本文算法提取的骨架線和骨架點(diǎn)數(shù)量中等，交叉點(diǎn)較少，拓?fù)溥B接數(shù)量錯(cuò)誤較少，并且受到噪聲影響時(shí)，提取的骨架點(diǎn)數(shù)量和拓?fù)溥B接數(shù)量變化不大，表現(xiàn)出了對(duì)噪聲的魯棒性。

表1 幾種方法的骨架點(diǎn)數(shù)量、拓?fù)溥B接數(shù)量、交叉點(diǎn)數(shù)量、拓?fù)溥B接錯(cuò)誤數(shù)量的統(tǒng)計(jì)Tab.1 Statistics of the number of skeleton points，topological connections，junction points and topological connection errors of several methods

3.3 算法效率

利用相同的軟硬件設(shè)備運(yùn)行L1中軸骨架線提取算法、MdCS算法、Point2Skeleton算法和本文算法，并對(duì)計(jì)算耗時(shí)進(jìn)行統(tǒng)計(jì)，結(jié)果如表2所示。從表2中可以看到，本文算法的時(shí)間損耗小于其他三種算法，即：略小于Point2Skeleton算法，相較于L1中軸骨架算法速度提升了10倍以上，比MdCS算法速度提升了百倍以上。由于Point2Skeleton是基于深度學(xué)習(xí)的算法，所以在計(jì)算骨架時(shí)相較于傳統(tǒng)算法更有優(yōu)勢(shì)，而該方法實(shí)際耗時(shí)表現(xiàn)略差與本文方法的原因是，Point2Skeleton在加載深度學(xué)習(xí)模型和保存運(yùn)算結(jié)果步驟中耗時(shí)較多。

表2 各算法的運(yùn)行時(shí)間Tab.2 Running time of different algorithms（s）

本文算法速度較快的原因是，使用了降采樣的方式降低了參與運(yùn)算的點(diǎn)的數(shù)量。骨架點(diǎn)迭代時(shí)，由于骨架點(diǎn)平衡狀態(tài)判斷的同時(shí)，也會(huì)在常數(shù)時(shí)間內(nèi)獲取骨架點(diǎn)迭代的方向，這種設(shè)計(jì)加速了骨架點(diǎn)的平衡移動(dòng)過(guò)程。

3.4 應(yīng)用

骨架線常被應(yīng)用于動(dòng)作捕捉識(shí)別［15-17］，基于骨架線可以改變模型的姿態(tài)和動(dòng)作，對(duì)于點(diǎn)云重建、數(shù)據(jù)生成以及形狀變化具有重要作用。直接通過(guò)三維圖形編輯軟件制作動(dòng)畫(huà)的效率瓶頸是導(dǎo)入的模型需要手動(dòng)綁定系統(tǒng)內(nèi)預(yù)置的骨架，包括骨架的三維坐標(biāo)和拓?fù)溥B接關(guān)系。此過(guò)程不僅耗費(fèi)時(shí)間，并且因?yàn)橄到y(tǒng)預(yù)置骨架的種類(lèi)有限，也會(huì)影響軟件的通用性。通過(guò)骨架線提取算法得到骨架線后，三維圖形圖像編輯軟件可以自動(dòng)獲取骨架點(diǎn)的三維坐標(biāo)、構(gòu)建拓?fù)溥B接，不再需要手動(dòng)調(diào)整參數(shù)和綁定骨架設(shè)定，可以提高動(dòng)畫(huà)制作的效率。

圖9（a）中展示的是，使用FARO Focus Premium激光掃描儀得到玩具馬（右上黃棕色）的三維點(diǎn)云模型（右下）；通過(guò)本文的算法提取到計(jì)算掃描點(diǎn)云的骨架線，如圖9（b）所示；將骨架線和點(diǎn)云模型導(dǎo)入三維圖形編輯軟件Blender中，可以輕松實(shí)現(xiàn)對(duì)于模型動(dòng)作的改變，圖9（c）中展示的是對(duì)馬的激光掃描點(diǎn)云基于骨架線的動(dòng)作變形，將馬的右前腿從彎曲動(dòng)作改變?yōu)橄蚯奥N起。圖1（e）和1（f）展示了基于骨架線將人模型的右腳和右手臂抬起后的效果（彩圖見(jiàn)期刊電子版）。

圖9 基于本文骨架線指導(dǎo)的激光點(diǎn)云的動(dòng)作變換Fig.9 Action transformation of lidar point cloud based on the guidance of curve skeleton

4 結(jié)論

本文提出了一種適用于激光掃描點(diǎn)云的骨架線的方法，并在此骨架基礎(chǔ)上實(shí)現(xiàn)了對(duì)模型的變形。通過(guò)將骨架點(diǎn)定義為具有對(duì)稱(chēng)法向且距離最近相關(guān)點(diǎn)的中點(diǎn)，通過(guò)迭代使中點(diǎn)平衡化得到骨架點(diǎn)。然后采取主成分分析法尋找滿(mǎn)足方向一致性要求的曲線段，使用廣度優(yōu)先搜索算法合并各曲線段，最后將拉普拉斯平滑后的曲線段相連得到完整的骨架線。通過(guò)在多個(gè)點(diǎn)云模型上開(kāi)展對(duì)比實(shí)驗(yàn)，驗(yàn)證了本文算法的有效性、魯棒性和高效性，將原本處理由33 041個(gè)點(diǎn)組成的點(diǎn)云所需86.110 s縮短到4.356 s。將提取的骨架應(yīng)用在激光雷達(dá)點(diǎn)云的動(dòng)作變換任務(wù)中，展示骨架線的實(shí)際作用。