三維點云法向量的模糊估值算法*

劉艷菊 張永德 楊波

(1.哈爾濱理工大學(xué) 智能機(jī)械研究所，黑龍江 哈爾濱 150080;2.齊齊哈爾大學(xué) 計算中心，黑龍江 齊齊哈爾 161006)

應(yīng)用三維點云的曲面重建技術(shù)有很多，但大多數(shù)技術(shù)在三維重建過程中都要求點云數(shù)據(jù)含有法向量.法向量是點云數(shù)據(jù)的一個基本屬性，用于標(biāo)識曲面的內(nèi)/外側(cè).雖然法向量可以通過三維激光掃描儀獲得的深度圖像而得到，但掃描過程中不可避免的噪聲、振動等因素使得到的法向量不準(zhǔn)確甚至有缺失，很難滿足曲面重建的要求，有必要對法向量重新進(jìn)行估值.

目前法向量估值算法總體分為兩類:局部估值和全局估值.在局部估值算法中，文獻(xiàn)［1］的原則組元分析(PCA)法通過計算兩點間歐式距離得到法向量所在直線，但該方法只適用于規(guī)則采樣表面，要求點云必須均勻，實際上很難達(dá)到這樣的要求.文獻(xiàn)［2］在此方法基礎(chǔ)上提出局部權(quán)值優(yōu)化投影因子(WLOP)算法，從原始數(shù)據(jù)中篩選提取能夠滿足PCA 要求的點云數(shù)據(jù)，然后進(jìn)行法向量估值，雖然效果較好但程序復(fù)雜度增加一倍.文獻(xiàn)［3］用Mahlanobis 距離替代歐式距離的算法，可以同時求出法向量所在直線及其方向，但對點云數(shù)據(jù)要求也非常高，同樣難以實現(xiàn).全局估值法向量的算法主要是由輸入點云數(shù)據(jù)構(gòu)建Voronoi 圖［4］，將Voronoi 圖中的每個極點看做是數(shù)據(jù)點的法向量;這類算法可以簡單高效地估算出法向量，但如果點云數(shù)據(jù)中含有噪聲點估算出來的法向量就會有較大的偏差.文獻(xiàn)［5］將局部算法和全局算法相結(jié)合，具有二者的優(yōu)勢但效率較低.

近期又有一些研究新思路:文獻(xiàn)［6］利用魯棒統(tǒng)計算法進(jìn)行法向量估值，該算法適用于含有尖銳特征的模型，而且需要人工調(diào)整參數(shù).文獻(xiàn)［7］利用圖的優(yōu)化理論對法向量的方向進(jìn)行推導(dǎo)，但此算法不適用于薄片類模型.文獻(xiàn)［8］采用變分算法既可以估算法向量所在直線又可以確定其方向，該算法對較薄模型推導(dǎo)效果較好，但對于含有尖銳特征的曲面容易產(chǎn)生錯誤結(jié)果.總體上，每種算法都在某些模型上能夠?qū)崿F(xiàn)較好的法向量估值而忽略其他情況，不能適用于任意形狀模型的法向量估值.

由于不同模型差別很大，很難用一個統(tǒng)一的數(shù)學(xué)模型表達(dá)出來.筆者在文獻(xiàn)［9］的基礎(chǔ)上提出基于模糊推理的法向量估值算法，對獲取的點云數(shù)據(jù)進(jìn)行單元劃分，并且求出每個單元中數(shù)據(jù)點的k 鄰近距離和曲率，將這兩項數(shù)據(jù)輸入到模糊系統(tǒng)中，根據(jù)人類經(jīng)驗給出模糊規(guī)則，通過模糊推理確定該部分?jǐn)?shù)據(jù)點適合的估值算法，然后按照相應(yīng)的法向量估值算法進(jìn)行計算.

1 法向量模糊估值系統(tǒng)及相應(yīng)算法

1.1 模糊推理系統(tǒng)

模糊推理系統(tǒng)［9］是由IF-THEN 規(guī)則構(gòu)建而成的，可以通過系統(tǒng)的方法將專家的語言信息結(jié)合到系統(tǒng)中.該系統(tǒng)包括模糊化、模糊規(guī)則集合、模糊推理機(jī)制和反模糊化，如圖1 所示.模糊化是將非模糊的輸入空間映射到模糊集合;反模糊化是執(zhí)行相反映射從模糊集合映射到非模糊的輸出空間;模糊集合的特征是由隸屬函數(shù)μ(·)表示;模糊推理機(jī)制是執(zhí)行基于模糊規(guī)則的一個從輸入的模糊集合空間U∈Rn向輸出模糊集合空間V∈Rm的映射.設(shè)U =U1U2…Un，Ui∈R，i =1，2，…，n;V =V1V2…Vm，Vk∈R，k=1，2，…，m.下面詳細(xì)介紹模糊推理系統(tǒng)在法向量估值上的具體應(yīng)用.

圖1 模糊推理系統(tǒng)組成Fig.1 Composition of fuzzy inference system

對于點云集合中某個數(shù)據(jù)點pk，其k 個相鄰的數(shù)據(jù)點集合為數(shù)據(jù)點的k 鄰近距離和曲率的組合能表示出模型所在區(qū)域的特征，所以將這兩個變量作為模糊規(guī)則的條件給出一組IF_THEN規(guī)則:

其中:x 是模糊推理系統(tǒng)的輸入向量，x =［xknn，xcur］T，xknn是點pk的k 個鄰近點的歐式距離之和，，距離之和大說明點云分布比較松散，距離之和小說明點云分布密集;xcur是點pk處的曲率;y 是模糊推理系統(tǒng)的輸出向量，(i =1，2)，Bk(k =1，2，…，m)分別是模糊集合Ui和Vk的語言變量，由隸屬函數(shù)(xi)和(yk)(j=1，2，…，M)表示.其中:xi是xknn或xcur中的任意一個，根據(jù)i 的取值來確定;yi是y 的某個取值;M 是模糊規(guī)則的總數(shù)量.使用推理機(jī)制、單值模糊化和中心平均反模糊化策略的模糊推理系統(tǒng)的輸出表示為

ξ(x)是模糊基函數(shù)向量，ξ(x)=［ξ1(x)，ξ2(x)，…，ξM(x)］T，其中ξj(x)可以被定義為

為了系統(tǒng)的計算效率，xknn和xcur的隸屬函數(shù)采用三角型隸屬函數(shù).

1.2 法向量的模糊估計算法

經(jīng)過對物體形狀的分析，總體上將物體形狀分成3類:局部平滑部分、薄片突出部分和不規(guī)則尖銳特征部分.一個物體可能含有其中1 種情況，也可能具有上述3 種情況.只需要將物體的點云數(shù)據(jù)中對應(yīng)的k 鄰近距離及其曲率輸入到模糊系統(tǒng)中，按照模糊規(guī)則就可以將此物體每一部分屬于哪種類型區(qū)分出來，根據(jù)不同類型給出不同的法向量估值算法.

對于平滑部分的法向量估值可以直接應(yīng)用PCA法計算法向量;對于薄片部分文中增加了一個檢測器，專門對此部分的法向量進(jìn)行特定估值;不規(guī)則尖銳部分采用增加附加點的方式進(jìn)行估值，具體方式為

其中，［1 α γ］是由模糊推理系統(tǒng)的輸出y 確定的.公式(4)中第1 項NPCA是每個數(shù)據(jù)點都要進(jìn)行的估算，所以其系數(shù)是1，對于平滑曲面只采用此估算就可以了.NPCA采用常用的原則組件分析方法，參見文獻(xiàn)［10］.Nchk和Natt是專門分別處理薄片部分和尖銳特征部分的算法.

1.2.1 檢測器估值算法

在曲率變化很大，而k 鄰近距離也很大時，說明該區(qū)域是屬于薄片特征的部分，可以增加檢測器［11］對法向量進(jìn)一步估值.檢測器沿兩個相鄰數(shù)據(jù)點的局部切線方向，按照k 鄰近搜索方式簡單地逼近薄片關(guān)鍵點.算法步驟如下:

(1)通過數(shù)據(jù)點pi的k 鄰近距離投影到其切平面上，來確定當(dāng)前pi的法向量的方向.

(2)如果pi的k 鄰近距離投影是在凸出部分的外側(cè)，那么pi被認(rèn)為是在薄片位置的數(shù)據(jù)點.

(3)取pi的k 鄰近點pj，(ni，nj)分別是這兩點的法向量所在直線，按照公式(5)推導(dǎo):

式中:Dij是pi、pj兩點的投影距離;fbs是線段的中點是沿著各自法向量所在直線距離pi和pj的單位長度;obs是fbs投影到線段的點.

(5)在薄片位置的數(shù)據(jù)點pj被重新定向，即ninj=-1，但它的方向不允許傳播給其他點.

(6)按照這種方式，在法向量方向迭代檢測過程中，執(zhí)行檢測器可以將相鄰的曲面方向翻轉(zhuǎn).

1.2.2 附加點估值算法

在曲率變化大且不規(guī)則、k 鄰近距離之和小的區(qū)域，點云數(shù)據(jù)分布密集，且含有尖銳特征.文中采用在曲面的內(nèi)側(cè)添加附加點［12］的方式幫助識別曲面的方向.附加點是從給定的點P 沿著法向量N 的方向投影產(chǎn)生的點，附加點與原始點距離是β，附加點集合={pn+1，pn+2，…，p2n}，按如下公式產(chǎn)生:

其中，N 是由PCA 方法計算每個數(shù)據(jù)點得到的無方向的法向量直線，N = {n1，n2，…，nn};通過設(shè)置β(0 ＜β≤0.001)的取值，附加點被定位在期望曲面的內(nèi)側(cè)，與原始數(shù)據(jù)非常接近;公式(6)中設(shè)pi是原始點云.將所有的數(shù)據(jù)點連接生成四面體模型，如圖2 所示.

圖2 加入附加點后的四面體Fig.2 Tetrahedral with attach points

在圖2 中，四面體可以分為3類，分別是:①四面體的各頂點只含有原始點云集合P 中的點;②四面體的各頂點同時含有集合P 和集合的點;③四面體的各頂點只含有集合中的點.在這3類四面體中，①和③類在法向量估計中是冗余的，可以刪掉.刪除①和③類四面體后，就能夠通過含有P 和的三角網(wǎng)格推導(dǎo)出恰當(dāng)?shù)姆较蛄?具體算法步驟如下:

(1)按照式(6)，求出附加點pn+i，構(gòu)建附加點集合.

(2)應(yīng)用合并集合中的數(shù)據(jù)點構(gòu)建四面體模型，具體方法見參考文獻(xiàn)［12-13］.

(3)在點云集合中，設(shè)一個數(shù)據(jù)點pk為起始點，構(gòu)建pk周圍的所有三角片集合{}，g 是環(huán)繞在點pk周圍的三角片的數(shù)量;構(gòu)建以pk為頂點的各角集合{}，如果角度之和為360°，pk點的法向量按公式(7)估算:

(4)構(gòu)建pk周圍的所有三角片邊長集合}，其中是第i 個三角片的3 條邊中的最長邊.因為冗余三角片與其他三角片相比總是有較長的邊，所以對三角片集合中各最長邊集合按升序排序.排序后求的和，j 是使等式成立的值.同時認(rèn)為是冗余三角片，將其刪除.接下來按式(7)求pk點的法向量值，然后執(zhí)行步驟(6).如果不存在j，按步驟(5)處理.

(5)采用加權(quán)求平均的方式計算法向量:

式中:hik為兩個數(shù)據(jù)點之間的距離，w(h)為權(quán)值函數(shù)、w(h)=1-6J2+8J3-4J4，J=h/R，R 為點pk的支撐半徑;pk的k 鄰近點相應(yīng)的法向量集合是已經(jīng)按照上面步驟(3)、(4)估算出的法向量的值;是所包含的鄰近點數(shù)量.

(6)循環(huán)檢測此區(qū)域的數(shù)據(jù)點，即可估算出相應(yīng)的法向量.

2 算法驗證

牙齒模型是一個典型的不規(guī)則模型(如圖3 所示)，牙齒的咬合面部分含有不規(guī)則的尖銳特征，切齒、犬齒含有薄片特征，而牙冠和牙齦部分的曲面較平滑.

圖3 不規(guī)則模型—牙頜Fig.3 Irregular model—denture

2.1 牙齒模型的點云數(shù)據(jù)分析

文中得到牙頜模型的點云數(shù)據(jù)后，計算每個數(shù)據(jù)的k 鄰近距離之和(i =1，2，…，n)，按升序排序并添加到優(yōu)先級隊列，從而確定點云的疏密程度.計算數(shù)據(jù)點的曲率xcur，在點云集中的區(qū)域沿著曲率變化最小的方向調(diào)整法向量的方向可以避免最大的可能誤差以及在具體實施過程中的“死鎖”問題，例如，含有尖銳特性的牙尖或牙窩.從xknn和xcur的取值可以知道點云數(shù)據(jù)密度和尖銳特征的分布情況，如圖4 所示.

圖4 點云的分布和曲率Fig.4 Distribution and curvature changes of point clouds

2.2 牙齒模型對應(yīng)的模糊規(guī)則

經(jīng)過對牙頜模型和實體牙齒的點云數(shù)據(jù)測量分析，給出相應(yīng)的模糊規(guī)則:

模糊規(guī)則中xknn和xcur的隸屬函數(shù)取值如圖5所示.圖5(a)中，“大”、“小”對應(yīng)的取值范圍分別是(0.05，+∞)、(-∞，0.05］;圖5(b)中“負(fù)無窮大”、“正無窮大”對應(yīng)的取值范圍分別是(-∞，-25］，［25，+∞)，“小”的范圍是(-25，25).

圖5 三角型隸屬函數(shù)Fig.5 Triangle membership functions

2.3 牙齒模型的法向量估值

文中采用21 號標(biāo)準(zhǔn)牙齒和牙頜模型的數(shù)據(jù)，按照前面給出的模糊推理系統(tǒng)，分別由切齒、犬齒、臼齒和牙頜模型的點云數(shù)據(jù)計算出相應(yīng)的k 鄰近距離和曲率，并輸入到模糊系統(tǒng)中，按照模糊推理規(guī)則將點云數(shù)據(jù)進(jìn)行分類估值法向量，其中計算含有尖銳特征時增加附加點使用的參數(shù)β =0.000 75.同時采用文獻(xiàn)［2］的WLOP 算法對上述數(shù)據(jù)進(jìn)行法向量估算.表1 給出文中算法與WLOP 算法執(zhí)行時間的比較結(jié)果.可以看出，數(shù)據(jù)量越大估值時間越長，當(dāng)模型中出現(xiàn)多種特征時，文中算法明顯優(yōu)于WLOP 算法的法向量估值.

表1 文中算法與WLOP 算法的執(zhí)行時間比較Table 1 Running time comparison between WLOP algorithm and the proposed one

兩種方法估算的法向量在進(jìn)行隱式曲面擬合［14-15］時，設(shè)擬合精度δ =5 ×10-4，圖6 分別給出上述4 種模型利用兩種法向量估值進(jìn)行曲面擬合的時間比較.從圖中可以看出，由于WLOP 算法法向量估算精度低，所以曲面擬合時間相對較長，尤其是對于數(shù)據(jù)量大的牙頜模型.曲面擬合的結(jié)果如圖7所示.

圖6 兩種算法的曲面擬合時間比較Fig.6 Comparison of surface fitting time between two algorithms

圖7 切牙、犬齒、臼齒和牙頜模型擬合結(jié)果Fig.7 Models of cutting，canine teeth，molar and denture

3 結(jié)語

文中提出了可對任意形狀物體的點云數(shù)據(jù)進(jìn)行法向量模糊估值的算法，利用模糊推理系統(tǒng)將點云數(shù)據(jù)進(jìn)行分類，將物體可能含有薄片特征和尖銳特征的點云區(qū)域區(qū)分處理:對于平滑的區(qū)域直接應(yīng)用PCA 方法估值其法向量，而含有薄片特征的區(qū)域由檢測器對關(guān)鍵點的法向量進(jìn)行推導(dǎo)，不規(guī)則尖銳特征的區(qū)域由增添附加點的算法對法向量估值.這樣使檢測器和附加點算法只應(yīng)用在局部區(qū)域，減少了計算量，而且針對不同區(qū)域應(yīng)用不同算法的法向量估值準(zhǔn)確率和曲面擬合效率均較高.

［1］Hoppe H，Tony D R，Tom D.Surface reconstruction from unorganized points［C］∥Proceedings of the 19th Annual Conference on Computer Graphics and Interactive Techniques.New York:ACM SIGGRAPH Computer Graphics，1992:71-78.

［2］Huang H，Li D，Zhang H，et al.Consolidation of unorganized point clouds for surface reconstruction ［J］.ACM Transactions on Graphics，2009，28(5):1-7.

［3］Jaakko L，Matthias Z，Emmanuel T，et al.A meshless hierarchical representation for light transport ［J］.ACM Transactions on Graphics，2008，27(3):1-37.

［4］Amenta N，Bern M.Surface reconstruction by Voronoi filtering［J］.Discrete ＆ Computational Geometry，1998，22(3):481-504.

［5］Alliez P，David C S，Tong Y Y，et al.Voronoi-based variational reconstruction of unoriented point sets ［C］∥Proceedings of the Fifth Eurographics Symposium on Geometry Processing.Barcelona:Eurographics Association，2007:39-48.

［6］Li B，Schnabel R W，Klein R，et al.Robust normal estimation for point clouds with sharp features［J］.Computer ＆Graphics，2010，34(2):94-106.

［7］Chen Y L，Lai S H.An orientation inference framework for surface reconstruction from unorganized point clouds［J］.IEEE Transactions on Image Processing，2011，20(3):762-775.

［8］Wang J，Yang Z W，Chen Fa-lai.A variational model for normal computation of point clouds［J］.The Visual Computer，2012，28(2):163-174.

［9］Liu Y J，Zhang Y D.A novel normal estimation based on fuzzy inference for mass point clouds of denture in the three-dimensional reconstruction ［J］.ICIC Express Letters Part B:Applications，2012，3(4):931-938.

［10］常俊彥，達(dá)飛鵬，蔡亮.基于特征融合的三維人臉識別［J］.東南大學(xué)學(xué)報:自然科學(xué)版，2011，41(1):47-51.Chang Jun-yan，Da Fei-peng，Cai Liang.3D face recognition based on feature fusion ［J］.Journal of Southeast University:Natural Science Edition，2012，41(1):47-51.

［11］Yukie N，Yutaka O，Hiromasa S.Smoothing of partition of unity implicit surfaces for noise robust surface reconstruction［J］.Computer Graphics Forum，2009，28(5):1339-1348.

［12］Taku I.Surface reconstruction with an interactive modification of point normals［C］∥Proceedings of the 2nd International Conference on Computational Modeling of Objects Represented in Images.Buffalo:CESER Publications，2010:263-274.

［13］Taku I，Yoshifumi K.Surface reconstruction from 3D scattered points with normals using both delaunay tetrahedralization and implicit function［C］∥Proceedings of the 7th Asia Simulation Conference on System Simulation and Scientific Computing.Beijing:Institute of Electrical and Electronics Engineers，2010:780-785.

［14］Stelldinger P，Latecki J L，Siqueira M.Topological equivalence between a 3D object and the reconstruction of its digital image［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，2007，29(1):126-140.

［15］Lange C，Polthier K.Anisotropic smoothing of point sets［J］.Computer Aided Geometric Dessign，2005，22(7):680-692.