朱厚盛，鮑憲帥，朱春元，陳明勝

(大連海事大學 信息科學技術學院，遼寧 大連 116026)

0 引 言

作為一種形狀描述的工具, 中軸概念[1]是由Blum于1967年首次提出并應用于圖像領域, 發(fā)展并推廣應用于高維情形。中軸可以看成是平面域內(nèi)邊界所有最大內(nèi)切圓的圓心集合[2],如圖1所示。由于中軸的良好特性, 因此被廣泛應用于各個領域中,例如形狀表示及變換[3]、醫(yī)學影像[4]、網(wǎng)絡規(guī)劃[5]、形狀分割[6]等。

圖1 最大圓模型

目前計算形狀中軸的方法普遍存在兩個缺點,即高時間成本或低中軸質量。首先,目前大多數(shù)方法都使用并行度不高的CPU作為主處理器,生成中軸的時間成本較高, 因而可以使用GPU的并行計算來快速計算出形狀中軸。其次, 目前中軸的計算方法離散程度較高, 導致求解出的中軸精度不足。

考慮以上因素, 本文提出了一種基于雙法線跟蹤的形狀中軸并行提取算法, 該算法可以高效得到形狀的高質量中軸。該算法步驟如下：①將形狀的包圍盒分割成若干的網(wǎng)格單元, 以便于接下來的并行化處理。②根據(jù)精度要求, 將形狀的邊界離散成由若干采樣點連接成的多邊形。③估計邊界點和邊界邊的法線方向。④以迭代方式計算每個采樣點的相應中軸點。最后,通過樣本點的拓撲連通性連接中軸點生成中軸。由于所有邊界點的法線方向都是高精度的,而使用這些法線來生成的中軸點是遵循中軸的定義的，因此所生成的中軸也是高質量的。

1 相關工作

按照對對象處理方式的不同, 可以將目前已存在的中軸提取算法劃大致分為三大類：中軸變換算法、細化算法和形狀分解算法。

中軸變換算法是Blum在1973年所提出的。中軸變換算法的主要原理是利用中軸的最大圓盤的定義, 即物體的中軸是物體內(nèi)部最大圓盤的圓心組成的點集, 由此可知中軸點與形狀邊界有兩個或兩個以上的切點。該算法主要用于處理多邊形模型或參數(shù)曲線曲面模型等連續(xù)模型對象, 通過利用嚴格的數(shù)學方法求解, 得到的結果精度較高, 計算過程較為復雜, 在處理復雜模型時穩(wěn)定性不高, 所以該方法常用于幾何建模和物體重建。Zhu等[7]提出了一種在平面域生成自由曲線邊界形狀中軸的方法。該方法在邊界曲線上采集足夠多的樣本點,利用這些樣本點計算Voronoi圖生成初始中軸；然后基于單向Hausdorff 距離對初始中軸進行裁剪, 得到平面域的最終中軸。Li[8]針對中軸變換的不穩(wěn)定性及冗余性提出了基于二次誤差最小化的中軸變換簡化算法(Q-MAT)。在此算法中將中軸變換的簡化轉換為中軸網(wǎng)格合并簡化的問題, 并對二次誤差度量框架進行了改進, 能夠有效的移除中軸變換中不重要的分支, 得到一個簡單緊湊且較為準確的中軸變換。而本文則同樣是利用網(wǎng)格來求取中軸, 主要是利用中軸定義跟蹤法線計算中軸, 并且加入了GPU并行計算, 從而大幅度提高了算法上的高效性。

細化算法一般是以像素為對象進行操作的。在數(shù)字圖像處理中, 我們通常認為當前像素點(物體圖像上的點)為黑色, 而背景像素點為白色。細化算法的主要原理就是利用迭代算法對圖像中的邊緣像素點進行剝離, 在保證拓撲性和連續(xù)性不變的情況下得到只有一個像素寬度的原圖像的簡化圖形的表示。此類算法應用于離散化對象, 如像素形狀、體素形狀、點云形狀等。Brunner等[9]提出了一種獲取形狀中軸的方法, 該方法是從模型邊界開始刪除無用的像素點, 并向內(nèi)部進行刪除操作, 直到?jīng)]有要刪除的點為止, 在這個過程中需要一些額外的因素來保持模型的拓撲特征。Hu等[10]在完全并行細化算法的基礎上提出了一種面向多孔介質的孔隙網(wǎng)絡中軸提取方法, 該方法計算較為簡單, 并且可以較好保持模型的全局特征, 但是該方法的重構過程較為復雜。Zhu等[11,12]提出了一種基于邊界擴張法的形狀中軸生成算法。由于離散方法通常是在離散對象上進行的,其精度低于中軸變換算法。

形狀分解算法的主要原理是先將原平面圖形進行分解, 再對分解后的小圖形求取中軸, 最后將得到的中軸連接起來合成原平面圖形的中軸。Wang等[13]提出了一種新的增量聚類算法, 此方法應用于書法字的中軸提取。首先通過K-均值聚類算法對書法字圖像進行分割,利用新的增量聚類算法對部分中軸進行提取,最后通過拓撲連接形成完整的書法字中軸。Zhu等[14]提出了一種利用擴張法將形狀平均劃分為不同部分并生成等分子中軸的方法, 考慮了相鄰部件之間的影響, 并結合這些子中軸點生成最終中軸。并行化可以提高算法的效率, 但由于被除數(shù)的限制, 并行度不高。

2 方法簡介

本文為了獲得形狀中軸, 首先將形狀的包圍盒分割成為若干的網(wǎng)格單元, 然后將形狀的邊界離散為若干樣本點連接成的多邊形。為了并行計算中軸點, 分別估計樣本點和邊界邊的內(nèi)法線方向(后文中所提及的法線皆為內(nèi)法線)。如圖2所示, 為了計算樣本點Pi的中軸點 (i=1,2,…,5), 需計算它們的內(nèi)切圓。根據(jù)中軸的定義可知, 內(nèi)切圓的圓心即是樣本點Pi的中軸點Mi。這里,內(nèi)切圓上的另一個邊界點Qi被稱為樣本點Pi的對偶邊界點,對偶邊界點Qi所在的邊界邊被稱為樣本點Pi的對偶邊界邊,樣本點Pi的內(nèi)切圓半徑ri被稱為樣本點的中軸半徑。對于任意一個樣本點P,為了計算其內(nèi)切圓, 需要首先搜索所有的邊界邊, 以便確定樣本點的對偶邊界邊。樣本點P的內(nèi)切圓可能和邊界有多個切點, 從而樣本點可以具有多于一個的對偶邊界邊, 但它們中的任何一個均足以計算出樣本點P的內(nèi)切圓。計算出對偶邊界邊后, 通過利用樣本點P及其對偶邊界邊的信息求解方程(詳情見5.3節(jié)), 可以計算出樣本點P的中軸點M和中軸半徑r。最后, 根據(jù)樣本點的拓撲關系, 連接所有樣本點的中軸點, 即可生成形狀中軸。

圖2 樣本點P1-P5的中軸點M1-M5和對偶邊界點Q1-Q5

本算法主要分成5個部分, 如圖3所示。

圖3 方法概述

(1)形狀分割：將形狀的包圍盒劃分成若干的網(wǎng)格單元, 每個網(wǎng)格單元的邊長為單元尺寸, 長度為gl。

(2)形狀離散：將形狀的邊界離散為若干樣本點連接成的多邊形。

(3)法線方向估計：估計邊界上所有樣本點及邊界邊的法線方向。

(4)中軸點計算：將所有的樣本點先放入待計算的樣本點集合中。在每一個追蹤回合i, 按照步驟(5)和步驟(6)計算中軸半徑在(i×gl,(i+1)×gl)的樣本點的中軸點, 并將已經(jīng)計算出中軸點的樣本點從集合中刪除。每個回合之后, 待計算的樣本點集合中的樣本點數(shù)量都會減少, 直到計算出所有樣本點的中軸點, 則中軸點計算過程結束。這里, 范圍(i×gl,(i+1)×gl)稱作回合i的當前范圍。

(5)樣本點的第一次法線跟蹤：對于每個未計算出中軸點的樣本點, 跟蹤它在當前范圍的法線段所在的網(wǎng)格單元, 將該樣本點存儲到網(wǎng)格中, 記錄為該網(wǎng)格的候選樣本點。每個樣本點的跟蹤需分配一個GPU線程。

(6)邊界邊的第二次法線跟蹤：對于每個邊界邊, 跟蹤它在當前范圍的法線段簇所在的網(wǎng)格單元。對于網(wǎng)格中每對候選樣本點和邊界邊, 計算出樣本點的中軸點。每個邊界邊的跟蹤需分配一個GPU線程。

(7)中軸生成：通過樣本點的拓撲連接性, 將每個樣本點對應中軸點連接起來形成中軸段, 進而形成形狀的中軸線。

3 前期工作

3.1 形狀分割

為了接下來的并行計算, 形狀所在的包圍盒需要進行分割成若干網(wǎng)格。通過切割x軸和y軸, 形狀的包圍盒被分割成多個正方形的網(wǎng)格單元, 其中網(wǎng)格單元的邊長被稱作單位長度gl。網(wǎng)格單元的數(shù)量太少會影響并行度, 太多則會增加跟蹤的迭代次數(shù), 因而可根據(jù)圖形大小適當分割。如圖4所示, 模型M在x軸和y軸上分別被分割6塊、7塊, 最終被分成42個網(wǎng)格單元。

圖4 形狀分割

將形狀分割后, 在后續(xù)的算法中, 可以將樣本點的搜索范圍限制于相應的網(wǎng)格單元中, 通過這種方法可以減少搜索范圍, 加快算法速度。同時,也便于引入GPU線程, 對不同網(wǎng)格的中軸點計算進行并行化處理。

3.2 形狀離散

在將形狀的包圍盒分割完成后, 需要對形狀的邊界進行離散。首先將形狀的邊界提取出來, 而后將邊界離散為若干等距離的樣本點, 并將樣本點相連接在一起。經(jīng)過離散操作后, 形狀就會離散成由這些樣本點連接而成的多邊形?？梢酝ㄟ^確定離散出的樣本點的數(shù)量來控制該方法的準確性和復雜性。樣本點數(shù)量越多, 最終得出的中軸越準確, 但相應的復雜性也會略有提高。

3.3 法向估計

根據(jù)中軸的定義, 樣本點的內(nèi)切圓與邊界至少有兩個切點。其中一個為樣本點本身, 另一個在某條邊界邊上。內(nèi)切圓的圓心為樣本點的中軸點, 中軸點就在兩個切點的法線上。為了計算每個樣本點的中軸點, 需要獲得樣本點和各個邊界邊的法向。

樣本點的法向是由和其相鄰的邊界邊的法向決定的。如圖5所示, 樣本點P點的相鄰邊界邊是AP與BP, 它們的法向NAP與NBP是分別垂直于邊界邊AP、BP的, 而P的法向NP則是NAP和NBP的中分線。

圖5 樣本點P的法向估計

邊界邊上的點的法向是由邊界邊端點的法向與該點在邊界邊的位置所共同決定的。如圖6所示,Q是邊界邊AB上的點,NA與NB分別是樣本點A、B的法向,設AQ∶QB=x∶y, 其中x+y=1, 則Q的坐標為：Q=xB+yA,Q的法向為

NQ=x×NB+y×NA

(1)

根據(jù)等式(1)可知, 當Q位于點A或B時,NQ等于NA或NB。當Q′為點Q附近一點時,Q′坐標為：Q′=(x+εx)B+(y+εy)A, 其中εx與εy趨近于0, 此時Q′的法向為

NQ′=(x+εx)×NB+(y+εy)×NA

(2)

比較等式(1)和等式(2)可知,Q′法向是無限接近于Q點法向的, 所以可知邊界邊上的點的法線方向是連續(xù)的。邊界邊上的點的法線的集合可以看成是該邊界邊的法線簇。圖6中,AS1和BS2長度為i×gl,S1E1和S2E2長度為gl, 那么線段E1S1和E2S2可以看成是點A和B在當前范圍的法線段, 而E1S1S2E2可以看成是AB在當前范圍的法線段簇。需要說明的是,在實際情況中,邊界邊AB的長度是較小的, 網(wǎng)格單元gl長度是較大的。文中這樣畫圖是為了方便說明、理解。

圖6 邊界邊上的點Q的法向估計

4 算法闡述

算法的主要步驟可以分為5個小步驟, 如圖7所示。

圖7 算法流程

(1)更新跟蹤的回合數(shù)i：在此步驟中對控制迭代的變量i進行加一操作。

(2)第一次法線跟蹤：跟蹤每一個未計算出中軸點的樣本點, 將樣本點的中軸點確定在某些網(wǎng)格單元中, 將樣本點記錄到網(wǎng)格單元中, 詳細步驟在5.1中具體描述。

(3)第二次法線跟蹤：跟蹤每一個邊界邊在當前范圍的法線段簇所在的網(wǎng)格單元, 在網(wǎng)格單元中記錄相應的邊界邊, 詳細步驟在5.2中具體描述。

(4)中軸點計算：在兩次跟蹤之后, 網(wǎng)格單元中會包含一些樣本點和邊界邊, 根據(jù)對應關系可計算出中軸點在該網(wǎng)格單元中的樣本點的中軸點, 詳細計算步驟在5.3中具體描述。

(5)判斷迭代結束條件：將已經(jīng)計算出中軸點的樣本點從樣本點集合中刪除, 如此時樣本點集合為空, 則算法迭代完成；如集合不為空, 則算法繼續(xù)。

5 雙法線跟蹤算法

根據(jù)中軸定義, 本文利用切點法向提出了一種并行雙法線跟蹤算法。通過該算法,可以迭代地計算出每個樣本點的中軸點, 進而得到形狀中軸。

雙法線跟蹤算法是一個迭代過程,涉及多回合跟蹤,其中第一次稱為第0回合。在第i回合跟蹤中,中軸半徑在(i×gl,(i+1)×gl)范圍內(nèi),即在當前范圍內(nèi)的樣本點將得到它們的中軸點。在每回合迭代之后,因為已經(jīng)能夠計算出一些樣本點對應的中軸點,所以這些樣本點不需要參加下一回合的計算,因而待計算的樣本點數(shù)量變小。經(jīng)過多回合迭代之后,所有的樣本點都得到了它們的中軸點,算法結束。在每一回合計算中,需要兩次法線跟蹤,即樣本點的第一次法線跟蹤和邊界邊的第二次法線跟蹤。

5.1 樣本點的第一法線跟蹤

樣本點的第一法線跟蹤是通過迭代地跟蹤其法向射線來計算它的中軸點。根據(jù)中軸定義可知，每個樣本點P的中軸點M必定在其法向射線NP上。如圖8所示，在第i回合中，設樣本點P的中軸半徑在當前范圍內(nèi)，則其中軸點M必定在當前范圍的法線段SE上。該法線段起始于點S，其中PS是i×gl，并在點E處結束，PE為(i+1)×gl。

圖8 樣本點P的法線追蹤

為了減少中軸點M的法線段的跟蹤區(qū)域，使用網(wǎng)格單元來跟蹤樣本點P的法線段。因為法線段SE包含中軸點M，所以只有與法線段SE相交的網(wǎng)格單元才可能包含中軸點。法線段SE的長度是gl，其起始點S和結束點E必定在同一網(wǎng)格單元或者邊相鄰或點相鄰的網(wǎng)格單元中。對于后兩種情況，最多有兩個或3個網(wǎng)格單元將分別與法線段SE相交,如圖9所示。這些與SE相交的網(wǎng)格單元將樣本點P記錄為其候選樣本點,以用于未來的中軸點計算。這就意味著，只有具有候選樣本點的網(wǎng)格單元才可能包含樣本點P的中軸點M。我們用鏈表來存儲網(wǎng)格單元的候選樣本點，以此最小化網(wǎng)格單元的存儲成本。

圖9 與法線段相交的網(wǎng)絡單元格

由于每個樣本點的第一次法線跟蹤沒有相互依賴性，可以獨立進行，所以可以充分利用GPU的并行計算特性，為每個樣本點分配一個GPU線程來并行實現(xiàn)。

5.2 邊界邊的第二法線跟蹤

在第一回合樣本點P的第一次法線跟蹤之后，它的中軸點M被鎖定在一些網(wǎng)格單元中，這些網(wǎng)格單元將該樣本點記錄為其候選樣本點。為了計算樣本點P的中軸點M，還要再計算出樣本點P的對偶邊界邊。由于每個邊界邊都可能是樣本點P的對偶邊界邊，因而需要檢查所有的邊界邊。

假定樣本點P的對偶邊界邊是T，對偶邊界點為Q。由于Q在邊界邊T上，其中軸點M一定位于邊界邊T的法線簇中。假設P的中軸半徑在(i×gl，(i+1)×gl)中，那么在第i回合中，中軸點M必定位于當前范圍的法線段簇中。因而，只有與該法線段簇相交的網(wǎng)格單元才可能包含中軸點M。如圖10所示，邊界邊AB的法線段簇與4個網(wǎng)格單元相交，因而中軸點M一定在這4個網(wǎng)格中，稱邊界邊AB是這些網(wǎng)格的候選邊界邊。對于每個網(wǎng)格，其可能不含有候選采樣點，也可能含有一個或者多個候選采樣點，則只有這些采樣點有可能與邊界邊T構成中軸點P。

圖10 與法線段簇相交的網(wǎng)格單元

在第i回合的兩次法線跟蹤后，對于每個網(wǎng)格都含有一些候選樣本點和候選邊界邊。如果中軸點M在某個網(wǎng)格單元中，則可以通過該網(wǎng)格單元中的每對候選樣本點和候選邊界邊的信息對，來計算出該中軸點。這個詳細過程在5.3節(jié)中闡述。

由于不同邊界邊的第二次法線跟蹤之間沒有耦合關系,所以也可以通過為每個邊界邊分配一個GPU線程的方式來并行實現(xiàn)。

5.3 中軸點與中軸生成

在第i回合的第二次法線跟蹤之后，一些網(wǎng)格單元中會含有一些候選樣本點和候選邊界邊對。對于某個網(wǎng)格單元，它在樣本點的第一次法線跟蹤時，記錄一些候選樣本點P1，P2，…，Pn。同時，它也在邊界邊的第二次法線跟蹤時，記錄邊界邊T1，T2，…，Tm,作為它的候選邊界邊。對于這些樣本點P1，P2，…，Pn,它們可能與邊界邊T1，T2，…，Tm上的任意一個點生成位于該網(wǎng)格單元中的中軸點M。假設中軸點是由樣本點P和邊界邊T生成出的，則根據(jù)中軸定義可知，中軸點M是切點為樣本點P和邊界邊T上的某一點的內(nèi)切圓的圓心，所以下面的等式成立

M=P+NP×r=Q+NQ×r

(3)

這里，r是中軸半徑，而點Q是樣本點P在邊界邊T上的對偶邊界點。通過等式(1)代替NQ，可將等式(3)轉化成如下的等式

(4)

等式(4)可以轉換成帶有未知量x，y的如下等式

(5)

這里a1、b1、c1、d1是通過P、A和B的坐標而得出的常量。而x+y=1，此時可以通過等式(5)得到如下等式

(6)

由于等式(6)在x軸和y軸分別成立，因此可以轉換成一個含有未知量r和未知量x的含有兩個方程的二元方程組，其解可以方便而準確地解出。注意，r的解的數(shù)量可能多于一個，但只采納中軸半徑在當前范圍的r。如果x、y和r都是非負的，則說明樣本點P的對偶邊界點Q在邊界邊T上，也可以確定采樣點P的中軸點M。這里，一個樣本點P可能通過不同的解得出多于一個中軸點，但是只采納中軸半徑最小的中軸點。

每回合迭代后，都有部分樣本點計算出其中軸點，則下一回合無需這些采樣點參與。通過多回合的迭代，最終可以計算出所有樣本點的中軸點。根據(jù)樣本點的拓撲連通性，可以對相應中軸點進行連接。如樣本點A、B是相連的，則其對應中軸點M1和M2也需要相連。通過中軸點的連接，可以形成形狀的中軸線。

6 實驗與分析

為了驗證算法的有效性，采取了多種不同類型的形狀進行實驗。實驗的平臺是Visual Studio2010，電腦處理器是英特爾i7-8700，顯卡是英偉達GTX1050Ti。

實驗共分成8組：實驗1是計算由簡單直線構成的凹多邊形的中軸；實驗2是計算由圓弧和直線構成的規(guī)則多邊形的中軸；實驗3是計算由直線和圓弧構成的帶孔的規(guī)則多邊形的中軸；實驗4是計算由不規(guī)則的曲線構成的隨機不規(guī)則圖形的中軸；實驗5～實驗8則是計算較為復雜的圖形的中軸。實驗結果如圖11所示。

圖11 形狀中軸顯示

經(jīng)過多次的實驗，可以得知，對于不同類型的形狀來說，該方法都能準確找到其中軸，并有效地顯示出來，可以得知該方法是有效的、精確的。

將傳統(tǒng)的細化算法和形狀分解算法與本文的算法進行比較實驗，并將本文算法分別通過CPU和GPU進行比對實驗，實驗時間對比見表1。

表1 實驗數(shù)據(jù)

通過上面的表格,我們可以很明顯看到，本文的算法比其它的算法在運行時間上有很大的提升。并且在利用GPU進行并行計算時，由于兩次的法線跟蹤都是并行計算的，所以花費的時間比較少。雖然理論上GPU線程多達數(shù)百個，但由于CPU線程運行速度實際比GPU線程快得多，所以GPU實際提升的效率大概在10倍左右。由此可知本文的方法是高效的。

本算法適用范圍廣泛，不僅可以計算規(guī)則形狀的中軸，也可以計算不規(guī)則形狀的中軸。在樣本點較多的情況下，其離散后的多邊形與原圖形相似，因而其生成的中軸與真實中軸非常接近。對于只有直線段的規(guī)則圖形，離散后的多邊形與原始圖形完全一致，邊界邊的法線簇和樣本點的法線簇也與原始圖形完全一致。由于本文的中軸點是嚴格按照中軸定義生成的，因而這種情況下生成的中軸即是精確中軸。由此可知本文生成的中軸是高質量的。

7 結束語

在本文中，提出了一種利用雙法線跟蹤來并行計算形狀中軸的算法。首先，將模型分割成網(wǎng)格模型，并將形狀邊界進行離散化，估計樣本點及邊界邊的法向。然后，對于模型的每個樣本點，使用并行雙法線跟蹤算法生成其對應中軸點。最后，根據(jù)其對應樣本點的連接性連接中軸點來生成形狀中軸。本文的主要貢獻在于：①利用GPU來對樣本點的第一次法線跟蹤和邊界邊的第二次法線跟蹤進行并行計算，計算效率顯著提高。②通過將模型的包裝盒劃分成網(wǎng)格單元，來減少并行計算中每個網(wǎng)格單元中候選樣本點和候選邊界邊對的平均數(shù)量，降低了算法的時間復雜度。③將樣本點按照其中軸半徑來分組，在每回合迭代中，僅中軸半徑在當前范圍內(nèi)的樣本點參與中軸點計算，減少了并行計算中每個樣本點的相交網(wǎng)格單元的平均數(shù)量，提高了算法效率。

未來的工作將會集中在將該方法合理地拓展到三維空間中，并進一步研究如何高效地將中軸線生成中軸面。