趙絲喆, 王寬全, 袁永峰

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

?

基于GPU的勢(shì)能場(chǎng)骨架提取并行算法

趙絲喆, 王寬全, 袁永峰

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

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

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

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

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

1勢(shì)能場(chǎng)骨架提取方法

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

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

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

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

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

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

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

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

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

2基于GPU的勢(shì)能場(chǎng)骨架提取并行算法

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

2.1并行性分析

(1)

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

2.2算法流程

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

圖1　GPU并行的勢(shì)能場(chǎng)骨架提取流程圖

2.3算法改進(jìn)

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

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

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

2.4線程分配

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

(2)

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

由此可知,對(duì)于不同的顯卡設(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ā)平臺(tái)為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ù)來(lái)自于文獻(xiàn)[6],分別是含有91 329個(gè)邊界點(diǎn),大小為204×132×260的結(jié)腸模型;含有6 555個(gè)邊界點(diǎn),大小為85×31×54的牛模型;含有9311個(gè)邊界點(diǎn),大小為87×74×45的螺旋模型;含有6 986個(gè)邊界點(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)行測(cè)試,均得到了較好的加速效果,說明本算法不受顯卡設(shè)備的制約．算法的加速比隨著數(shù)據(jù)規(guī)模的增大而逐漸提升,當(dāng)數(shù)據(jù)規(guī)模達(dá)到256×256×487時(shí),普通并行算法在兩種顯卡上分別可以達(dá)到11倍和15倍的加速比,而加入了常量存儲(chǔ)器和共享存儲(chǔ)器的改進(jìn)算法又進(jìn)一步將加速比提升了20%左右,分別達(dá)到了15倍和18倍．

圖2　GPU并行的勢(shì)能場(chǎng)骨架提取結(jié)果

體數(shù)據(jù)規(guī)模邊界點(diǎn)個(gè)數(shù)串行算法運(yùn)行時(shí)間/s并行加速比(FX5800)簡(jiǎn)單并行加常量存儲(chǔ)器加共享存儲(chǔ)器并行加速比(GTX580)簡(jiǎn)單并行加常量存儲(chǔ)器加共享存儲(chǔ)器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線程選擇

分別對(duì)基于GT200(計(jì)算核心1.3)和GF100(計(jì)算核心2.0)架構(gòu)的FX5800和GTX580進(jìn)行了測(cè)試．由式(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)表明,對(duì)于FX5800,GPU的占用率對(duì)加速比產(chǎn)生了較大影響．當(dāng)Threads per Block為193時(shí),GPU占用率從93.8%下降到87.5%,而加速比也相應(yīng)從14.1下降到了10.5．隨后當(dāng)Threads per Block的值為321時(shí),又再次出現(xiàn)了加速比隨GPU占用率的降低而急速下降的情況．而圖3(b)表示,對(duì)于GTX580而言,這種關(guān)聯(lián)并不明顯,這是由于GF100架構(gòu)中新增加了高速緩存功能,降低了訪問沖突的發(fā)生．但仍可以發(fā)現(xiàn),幾個(gè)加速比的低峰值,均出現(xiàn)在GPU占用率較小的情況下．因此,在分配線程時(shí)應(yīng)選擇合適的Threads per Block,使得GPU的占用率足夠大,從而得到更好的加速比．

(a) FX5800

(b) GTX580

4結(jié)語(yǔ)

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

參考文獻(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

中圖分類號(hào):P315.69

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

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

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

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

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

收稿日期:2015-02-03.

doi：10.11918/j.issn.0367-6234.2016.05.002