基于GPU的勢能場骨架提取并行算法

趙絲喆, 王寬全, 袁永峰

(哈爾濱工業(yè)大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院, 150001 哈爾濱)

?

基于GPU的勢能場骨架提取并行算法

趙絲喆, 王寬全, 袁永峰

(哈爾濱工業(yè)大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院, 150001 哈爾濱)

摘要:為解決勢能場骨架提取方法計(jì)算效率低、提取過程耗時大的問題,同時為降低該方法的時間復(fù)雜度,提出了基于GPU的勢能場骨架提取并行算法,并充分利用CUDA架構(gòu)特有的常量存儲器和共享存儲器對普通并行算法進(jìn)行改進(jìn)．討論了如何根據(jù)程序和顯卡設(shè)備的固有屬性來分配線程以達(dá)到最高的GPU占用率,從而得到最優(yōu)的加速效果．對多組3D模型進(jìn)行測試的結(jié)果表明,隨著數(shù)據(jù)規(guī)模的增大,加速效果逐漸提升,處理256×256×487的體數(shù)據(jù)時,可獲得18倍的加速比．

關(guān)鍵詞:圖形處理器;并行計(jì)算;勢能場;骨架提取;通用并行計(jì)算架構(gòu)

3D物體的骨架類似于二維情形下的中軸線,可以形象地表述物體的拓?fù)涮卣?被廣泛應(yīng)用于計(jì)算機(jī)動畫、機(jī)械設(shè)計(jì)、醫(yī)學(xué)圖像、虛擬導(dǎo)航和虛擬內(nèi)窺鏡等技術(shù)領(lǐng)域[1]．目前,骨架提取的方法可以根據(jù)輸入數(shù)據(jù)類型的不同來劃分．文獻(xiàn)[2]針對計(jì)算機(jī)動畫中的三角網(wǎng)格面數(shù)據(jù)計(jì)算出數(shù)據(jù)的voronoi圖,應(yīng)用平均曲率對其進(jìn)行邊緣劃分,進(jìn)而得到離散的骨架點(diǎn);文獻(xiàn)[3]針對三維激光掃描技術(shù)獲得的點(diǎn)云數(shù)據(jù)利用馬爾可夫隨機(jī)場模型的3D非剛性匹配技術(shù)追蹤點(diǎn)云數(shù)據(jù),得到的軌跡即為物體骨架;此外,還有CT、MRI等序列掃描圖像合成的體數(shù)據(jù),由于體數(shù)據(jù)可以完整保留物體的內(nèi)部信息,在醫(yī)學(xué)圖像可視化和生物信息等方面都具有重要意義．

針對于體數(shù)據(jù)的骨架提取方法主要包括:拓?fù)浼?xì)化法,如文獻(xiàn)[4]中反復(fù)剝離物體的最外層體素,直至單連通即為骨架;距離場方法,如文獻(xiàn)[5]對數(shù)據(jù)進(jìn)行距離變換,抽取出局部距離極大點(diǎn)作為骨架點(diǎn)．但文獻(xiàn)[6]表示以上兩種方法所提取出的骨架易受噪聲干擾．因此,Cornea ND等人提出了勢能場方法[7],綜合考慮到所有表面點(diǎn)的影響,而不像距離場只考慮距離最近的單個表面點(diǎn),因此該方法具備更好的魯棒性．但實(shí)驗(yàn)表明,運(yùn)用該方法作用于256×256×487的虛擬人體數(shù)據(jù)時,需要耗費(fèi)5 h,這在臨床應(yīng)用中是無法忍受的．因此,本文提出基于圖形處理器(GPU)的勢能場骨架提取并行算法,并加入CUDA架構(gòu)特有的存儲器結(jié)構(gòu)來優(yōu)化算法,討論了如何根據(jù)顯卡設(shè)備和程序的固有屬性來分配線程以達(dá)到最高的GPU占用率,從而得到更優(yōu)的加速效果．

1勢能場骨架提取方法

物體的骨架需要具備以下特性:1)纖細(xì)性,凝練出物體的基本結(jié)構(gòu);2)連通性,若物體本身是連通的,則骨架也同樣是連通的;3)中心性,骨架應(yīng)處于物體的核心位置上．針對體數(shù)據(jù)而言,骨架指的是由單體素組成的、連通的、位于物體中心的一序列體素[8]．

Cornea ND提出的勢能場方法假設(shè)物體表面遍布同種點(diǎn)電荷,作用于物體內(nèi)部,從而形成靜電斥力場．通過計(jì)算各點(diǎn)在勢能場中的受力情況,選取特殊的點(diǎn)作為種子點(diǎn),連接這些點(diǎn)形成3D物體的核心骨架．具體步驟如下:

1)將體數(shù)據(jù)中的各個離散體素點(diǎn)分類為外部點(diǎn)、表面點(diǎn)、邊界點(diǎn)和內(nèi)部點(diǎn)．外部點(diǎn)是體素值為0的點(diǎn);表面點(diǎn)是指其26鄰域中至少存在一個外部點(diǎn);邊界點(diǎn)是指其26鄰域中至少存在一個表面點(diǎn);其余的均為內(nèi)部點(diǎn)．

2)勢能場的計(jì)算．內(nèi)部點(diǎn)或邊界點(diǎn)P會受到周圍的表面點(diǎn)C所帶來的斥力,該斥力與距離成反比,計(jì)算公式為

3)選取場值為0且場方向發(fā)生改變的點(diǎn)為關(guān)鍵點(diǎn)．首先檢測一個立方體區(qū)域內(nèi)的8個頂點(diǎn),若場值在X,Y,Z方向上均發(fā)生改變,則該區(qū)域中可能存在關(guān)鍵點(diǎn),迭代劃分立方體,直至得到一個不可分的體素點(diǎn)為止．接著計(jì)算該點(diǎn)勢場力的雅克比矩陣,根據(jù)以下定義將關(guān)鍵點(diǎn)分類:

定義1雅克比矩陣特征值的實(shí)部和虛部均為負(fù)時,該點(diǎn)稱為吸引點(diǎn),其周圍所有向量都指向該點(diǎn)．

定義2雅克比矩陣特征值的實(shí)部和虛部均為正時,該點(diǎn)稱為排斥點(diǎn),其周圍所有向量都背離該點(diǎn)．

定義3雅克比矩陣特征值的實(shí)部和虛部有正有負(fù)時,該點(diǎn)稱為鞍點(diǎn),其周圍向量有的指向該點(diǎn),有的背離該點(diǎn)．

4)骨架生長和骨架連接．遍歷所有鞍點(diǎn),每個鞍點(diǎn)均生成一個骨架段．從鞍點(diǎn)出發(fā),以正特征值對應(yīng)的特征向量為方向,使用力跟隨法按照一定的步長前進(jìn)．連接所有鞍點(diǎn)形成骨架．

2基于GPU的勢能場骨架提取并行算法

通過實(shí)際運(yùn)行發(fā)現(xiàn),勢能場的計(jì)算占據(jù)了整個提取過程98%的時間,因此,減少骨架提取時間的核心思想就是加快勢能場的計(jì)算．

2.1并行性分析

(1)

其中,i為內(nèi)部點(diǎn),j為邊界點(diǎn),k為表面點(diǎn),勢能場計(jì)算的時間復(fù)雜度是Ο((i+j)×k)．由式(1)可知,勢場力的計(jì)算具備獨(dú)立性,點(diǎn)與點(diǎn)之間不互相影響,因此可以將原本的串行計(jì)算并行化．通常情況下,內(nèi)部點(diǎn)和邊界點(diǎn)的個數(shù)大于表面點(diǎn)的個數(shù),外層循環(huán)的計(jì)算量更大,于是本算法將外層循環(huán)放至GPU中,為每一個內(nèi)部點(diǎn)和邊界點(diǎn)分配一個線程,將時間復(fù)雜度降至Ο(k),且k遠(yuǎn)小于(i+j)．另外,根據(jù)定義可知邊界點(diǎn)處于表面點(diǎn)和內(nèi)部點(diǎn)之間,而計(jì)算平均值要比計(jì)算距離簡單很多,于是本算法將邊界點(diǎn)26鄰域點(diǎn)的平均勢場力作為該點(diǎn)的勢場力,此過程同樣放至GPU中．

2.2算法流程

在基于GPU的并行系統(tǒng)中,CPU和GPU各司其職,CPU負(fù)責(zé)復(fù)雜的流控制等需要串行處理的部分,而密集型數(shù)據(jù)的并行計(jì)算部分則交由GPU完成[9]．有別于原始的串行算法,本文中CPU只負(fù)責(zé)簡單的體素點(diǎn)劃分,而復(fù)雜的勢場力計(jì)算則交由GPU完成．計(jì)算時將各個內(nèi)部點(diǎn)和邊界點(diǎn)平均分配給每個線程,多線程并行執(zhí)行．具體的程序流程見圖1．

圖1　GPU并行的勢能場骨架提取流程圖

2.3算法改進(jìn)

在CUDA架構(gòu)中存在多種存儲結(jié)構(gòu),按照存取速度由快到慢排列依次是:寄存器(Register)、常量存儲器(Constant Memory)、共享存儲器(Shared Memory)、紋理存儲器(Texture Memory)、局部存儲器(Local Memory)和全局存儲器(Global Memory)．

在普通的并行算法中,僅僅使用了寄存器和全局存儲器,為了進(jìn)一步加快程序的運(yùn)行速度,提出了改進(jìn)的并行算法．結(jié)合算法本身的特點(diǎn),使用了訪問速度可以與寄存器媲美的常量存儲器和共享存儲器．

首先,對于每個內(nèi)部點(diǎn)而言,周圍的表面點(diǎn)個數(shù)有限,數(shù)據(jù)量小于常量存儲器的容量64 kB,并且表面點(diǎn)僅用于讀取,并不對其改寫,因此可以將表面點(diǎn)的信息存儲在常量存儲器中以加快內(nèi)部點(diǎn)的勢場力計(jì)算．另外,在程序執(zhí)行過程中,需要頻繁地訪問全局存儲器來修改勢場力數(shù)組的值,這無疑會帶來較大的延遲,因此在改進(jìn)算法中將勢場力數(shù)組移入到共享存儲器中．

2.4線程分配

對于串行程序而言,程序的執(zhí)行時間可以表達(dá)為關(guān)于問題規(guī)模(即輸入數(shù)據(jù)規(guī)模)的函數(shù),而對于并行程序而言,執(zhí)行時間不僅與問題規(guī)模相關(guān),還與并行體系結(jié)構(gòu)相關(guān),而GPU的占用率(occupancy)就是其中一個重要的考慮參數(shù)[10],計(jì)算公式如下MP(Multiprocessor):

(2)

maxRegister per MP,maxSharedMemory per MP和maxThreadBlocks per MP分別表示最大寄存器數(shù)量、最大共享存儲器容量以及最大線程塊數(shù)目,其值均取決于顯卡的規(guī)格參數(shù);而Registers per Thread和SharedMemory per Block則表示程序運(yùn)行中實(shí)際用到的寄存器數(shù)量和共享存儲器的容量,其值均由CUDA在啟動線程時自動分配．

由此可知,對于不同的顯卡設(shè)備,可先將顯卡和程序本身的固有屬性參數(shù)代入式(2),計(jì)算出合適的Threads per Block值即每塊線程數(shù),以達(dá)到最優(yōu)的GPU占用率,從而得到更好的加速效果．

3實(shí)驗(yàn)結(jié)果與性能分析

3.1實(shí)驗(yàn)環(huán)境

開發(fā)平臺為3.33 GHz Intel(R) Xeon(R) X5680 雙核CPU,24G內(nèi)存,顯卡分別為NVIDIA Quadro FX 5800和NVIDIA GeForce GTX 580,顯卡的各項(xiàng)參數(shù)如表1所示．開發(fā)工具為Microsoft Visual Studio 2010和CUDA4.0．

表1　顯卡的規(guī)格參數(shù)(MP: Multiprocessor)

3.2提取結(jié)果

圖2顯示了本程序的骨架提取結(jié)果,采用的體數(shù)據(jù)來自于文獻(xiàn)[6],分別是含有91 329個邊界點(diǎn),大小為204×132×260的結(jié)腸模型;含有6 555個邊界點(diǎn),大小為85×31×54的牛模型;含有9311個邊界點(diǎn),大小為87×74×45的螺旋模型;含有6 986個邊界點(diǎn),大小為54×87×75的恐龍模型．

3.3加速比

在不同的輸入數(shù)據(jù)規(guī)模下,普通并行算法和改進(jìn)并行算法的加速比結(jié)果如表2所示．將Threads per Block和Blocks per Grid均設(shè)定為128,分別在顯卡FX5800和GTX580上進(jìn)行測試,均得到了較好的加速效果,說明本算法不受顯卡設(shè)備的制約．算法的加速比隨著數(shù)據(jù)規(guī)模的增大而逐漸提升,當(dāng)數(shù)據(jù)規(guī)模達(dá)到256×256×487時,普通并行算法在兩種顯卡上分別可以達(dá)到11倍和15倍的加速比,而加入了常量存儲器和共享存儲器的改進(jìn)算法又進(jìn)一步將加速比提升了20%左右,分別達(dá)到了15倍和18倍．

圖2　GPU并行的勢能場骨架提取結(jié)果

體數(shù)據(jù)規(guī)模邊界點(diǎn)個數(shù)串行算法運(yùn)行時間/s并行加速比(FX5800)簡單并行加常量存儲器加共享存儲器并行加速比(GTX580)簡單并行加常量存儲器加共享存儲器16×16×487109713.13.0x3.3x3.6x2.7x2.4x2.6x32×32×4872225424.85.3x5.8x6.2x5.0x5.3x5.6x64×64×48745464213.67.2x7.8x8.8x6.4x7.0x7.4x128×128×4871016191800.98.8x9.6x11.1x11.0x12.1x13.0x256×256×48720601016212.211.2x12.4x14.7x15.2x16.9x18.7x

3.4線程選擇

分別對基于GT200(計(jì)算核心1.3)和GF100(計(jì)算核心2.0)架構(gòu)的FX5800和GTX580進(jìn)行了測試．由式(2)可知,影響GPU占用率的主要參數(shù)是Registers per Thread,SharedMemory per Block和Threads per Block．運(yùn)行程序前,在編譯選項(xiàng)里面加入--ptxas-options=-v命令,可以查看到本程序的Registers per Thread的值為42,SharedMemory per Block的值為1548字節(jié)．將參數(shù)代入式(2),再調(diào)節(jié)Threads per Block的值,就可以得到不同的GPU占用率,如圖3所示．圖3(a)表明,對于FX5800,GPU的占用率對加速比產(chǎn)生了較大影響．當(dāng)Threads per Block為193時,GPU占用率從93.8%下降到87.5%,而加速比也相應(yīng)從14.1下降到了10.5．隨后當(dāng)Threads per Block的值為321時,又再次出現(xiàn)了加速比隨GPU占用率的降低而急速下降的情況．而圖3(b)表示,對于GTX580而言,這種關(guān)聯(lián)并不明顯,這是由于GF100架構(gòu)中新增加了高速緩存功能,降低了訪問沖突的發(fā)生．但仍可以發(fā)現(xiàn),幾個加速比的低峰值,均出現(xiàn)在GPU占用率較小的情況下．因此,在分配線程時應(yīng)選擇合適的Threads per Block,使得GPU的占用率足夠大,從而得到更好的加速比．

(a) FX5800

(b) GTX580

4結(jié)語

針對勢能場方法提取骨架耗時長、計(jì)算復(fù)雜度高等問題,提出了一種基于GPU的勢能場骨架提取并行算法．由于在提取過程中勢能場的計(jì)算占據(jù)了98%的時間,對勢能場的計(jì)算進(jìn)行了并行性分析并提出并行算法,另加入了CUDA架構(gòu)特有的存儲器結(jié)構(gòu)對普通并行算法進(jìn)行了優(yōu)化．此外,還討論了如何根據(jù)顯卡設(shè)備和程序的固有屬性來分配線程以達(dá)到最高的GPU占用率,從而得到最優(yōu)的加速效果．實(shí)驗(yàn)結(jié)果表明,當(dāng)處理256×256×487規(guī)模的體數(shù)據(jù)時,可以獲得18倍的加速比,有效地減少了骨架提取時間．

參考文獻(xiàn)

[1] LIVESU M, SCATENI R. Extracting curve-skeletons from digital shapes using occluding contours[J]. The Visual computer, 2013, 29(9):907-916.

[2] TAGLIASACCHI A, ALHASHIM I,OLSON M,et al. Mean Curvature Skeletons[J]. Computer Graphics Forum, 2012, 31:1735-1744.

[3] Zhang Q, SONG X,SHAO X,et al. Unsupervised skeleton extraction and motion capture from 3D deformable matching[J]. Neurocomputing, 2013,100:170-182.

[4] BERTRAND G, COUPRIE M. Powerful parallel and symmetric 3d thinning schemesbased on critical kernels[J]. J Math Imaging Vis, 2014, 48:134-148.

[5] ARCELLI C, DI BAJA G S, SERINO L. Distance-driven skeletonization in voxel images[J]. IEEE Trans Pattern Anal Mach Intell, 2011, 33(4):709-720.

[6] CORNEA ND, SILVER D, MIN P. Curve-skeleton properties, applications, and algorithms[J]. IEEE Transactions on Visualization and Computer Graphics,2007,13(3):530-548.

[7] CORNEA N D, SILVER D,YUAN Xiaosong,et al.Computing hierarchical curve-skeletons of 3D objects[J].The Visual Computer,2005,21(11):945-955.

[8] LIU B, TELEA A C,ROERDINK J B T M,et al. Parallel centerline extraction on the GPU[J]. Computers & Graphics, 2014,41:72-83.

[9] LEE C, RO W W,GAUDIOT J L,et al. Boosting CUDA applications with CPU-GPU hybrid computing[J]. International Journal of Parallel Programming, 2014, 42(2):384-404.

[10]JIMéNEZ J. Three-dimensional thinning algorithms on graphics processing units and multicore CPUs[J]. Concurrency and Computation: Practice and Experience, 2012, 24:1551-1571.

(編輯王小唯苗秀芝)

Parallel method of skeleton extraction using potential field on GPU

ZHAO Sizhe, WANG Kuanquan,YUAN Yongfeng

(School of Computer Science and Technology, Harbin Institute of Technology, 150001 Harbin,China)

Abstract:For curve skeleton extraction algorithm, in order to improve the efficiency of potential field computation and save the time of extraction process, we presented a parallel potential field skeleton extraction method to reduce the time complexity, which was suitable for implementation on GPU, and then improved it by using constant memory and shared memory which was unique in CUDA. In order to achieve the highest GPU occupancy and the best speedups, we discussed how to assign threads according to the property of program and graphics device. The implementation was tested on several complex 3D models in CUDA framework. The results showed that our method had excellent performance especially on large data scale. When processing the volume data with the scale of 256×256×487, this improved method achieved speedups of 18x.

Keywords:GPU; parallel computing; potential field; skeleton extraction; CUDA

中圖分類號:P315.69

文獻(xiàn)標(biāo)志碼:A

文章編號:0367-6234(2016)05-0018-05

通信作者:王寬全,wangkq@hit.edu.cn.

作者簡介:趙絲喆(1991—),女,碩士研究生;王寬全(1964—),男,教授,博士生導(dǎo)師.

基金項(xiàng)目:國家自然科學(xué)基金面上項(xiàng)目(61173086).

收稿日期:2015-02-03.

doi：10.11918/j.issn.0367-6234.2016.05.002