齊 金，李 寬，楊燦群，杜云飛

(1.國防科學(xué)技術(shù)大學(xué)并行與分布處理重點實驗室，湖南 長沙 410073;2.國防科學(xué)技術(shù)大學(xué)計算機學(xué)院，湖南 長沙 410073)

XeonPhi平臺上基于模板優(yōu)化的3DGVF場計算加速*

齊 金1，李 寬2，楊燦群1，杜云飛2

(1.國防科學(xué)技術(shù)大學(xué)并行與分布處理重點實驗室，湖南 長沙 410073;2.國防科學(xué)技術(shù)大學(xué)計算機學(xué)院，湖南 長沙 410073)

3D梯度向量流場(3D GVF field)廣泛應(yīng)用于多種3D圖像分析算法中，其計算需要多次迭代，計算量大，如何提高其計算速度具有重要的研究意義。面向Intel Xeon Phi眾核集成架構(gòu)，首次進(jìn)行了3D GVF場計算的加速優(yōu)化。首先，挖掘3D圖像像素點間存在的天然并行性，發(fā)揮眾核架構(gòu)優(yōu)勢，嘗試線程級并行(多核)和數(shù)據(jù)級并行(SIMD)。其次，3D GVF場的計算過程是一種典型的3D-7點模板運算，結(jié)合Xeon Phi架構(gòu)的L2 緩存規(guī)格，提出一種高效的數(shù)據(jù)分塊策略，充分挖掘數(shù)據(jù)的時/空局部性，有效緩解模板計算引起的緩存缺失，提升了計算性能。實驗結(jié)果表明，引入模板優(yōu)化技術(shù)能顯著提升3D GVF場的計算速度，在圖像維度為5123時，所提方法在57核Xeon Phi平臺上的性能相比在2.6 GHz 8核16線程的Intel Xeon E5-2670 CPU上的性能，加速比可達(dá)2.77。

3D梯度向量流場；Xeon Phi；模板優(yōu)化；緩存分塊

1 引言

在基于變分法的圖像處理中，主動輪廓模型(Active Contour Model)廣泛應(yīng)用于邊界檢測、圖像分割和運動跟蹤[1～3]。概括來說，主動輪廓模型將曲線(曲面)與能量函數(shù)相聯(lián)系，由內(nèi)力和外力場共同引導(dǎo)，使曲線(曲面)不斷向能量最小化的方向演化。其中，內(nèi)力由曲線決定，外力場由圖像計算得到。諸多的外力場定義中，Xu C等人[4]提出的梯度向量流場GVF場(Gradient Vector Flow Field)通過一組偏微分方程對圖像的梯度向量進(jìn)行擴散，具有捕獲范圍大、抗噪聲等優(yōu)點，成為主動輪廓模型中經(jīng)典的外力場，得到了廣泛的研究與應(yīng)用[5,6]。

GVF場的計算與圖像大小規(guī)模緊密相關(guān)，且需多次迭代。3D圖像數(shù)據(jù)量大，其3D GVF場的計算速度一直是制約其應(yīng)用的瓶頸，如何提高GVF場，尤其是3D GVF場的計算速度已獲得研究者的關(guān)注。在高性能計算領(lǐng)域，已有研究使用GPU加速3D GVF場計算，如GPU+OpenGL[7]、GPU+OpenCL[8]等。圖像中各像素的處理存在天然的并行性，加上GPU中獨特的紋理存儲(Texture Memory)能方便存取鄰居像素數(shù)據(jù)，3D GVF場在GPU平臺上取得了較好的加速效果。

Intel 新推出的Xeon Phi眾核集成架構(gòu)，提供大量IA(Intel Architecture)架構(gòu)的輕量核，在兼容傳統(tǒng)編程模型的基礎(chǔ)上，能提供更高的計算性能。在此背景下，本文首次面向Xeon Phi平臺進(jìn)行3D GVF場計算的并行優(yōu)化。主要的工作包括兩個方面：(1)挖掘3D圖像像素點間存在的天然并行性，發(fā)揮眾核架構(gòu)優(yōu)勢，嘗試線程級并行(多核)和數(shù)據(jù)級并行(SIMD)。并對圖像數(shù)據(jù)在Xeon Phi協(xié)處理器內(nèi)存中的存取模式進(jìn)行優(yōu)化，以節(jié)省指令、提高吞吐率，這對其他圖像算法在Xeon Phi上的處理均有借鑒意義。(2)計算3D GVF場時，需存取3D空間中的鄰居像素，這是一種典型的3D-7點Laplace模板計算。在Xeon Phi平臺上，結(jié)合L2緩存大小，提出一種高效、具體的分塊策略，充分挖掘數(shù)據(jù)的時/空局部性，有效緩解模板計算引起的緩存缺失，提升了計算性能。實驗結(jié)果表明,引入模板優(yōu)化技術(shù)能顯著提升3D GVF場的計算速度。

2 背景和相關(guān)工作

(1)3D GVF場計算方法。

3D GVF場是能使得如下能量函數(shù)最小的圖像空間的向量場V；

(1)

其中，V0是初始向量場。

上述最優(yōu)化問題可通過如下歐拉方程解得：

(2)

該問題對應(yīng)的數(shù)值解法如算法1所示：

算法13D GVF數(shù)值迭代解法

fori∈[1Iterations] do

forpoint(x,y,z)∈the image do

laplacian←-6Vi(x,y,z)+(x+1,y,z)+Vi(x-1,y,z)+(x,y+1,z)+Vi(x,y-1,z)+(x,y,z+1)+Vi(x,y,z-1)

Vi+1(x,y,z)←Vi(x,y,z)+μ×laplacian-(Vi(x,y,z)-V0(x,y,z))|V0(x,y,z)|2

end for

end for

可見，同次迭代中，各點的梯度向量更新是相對獨立的，這種天然的數(shù)據(jù)并行性能很好地發(fā)揮Xeon Phi的強大計算能力，為后續(xù)的線程級并行(多核)和數(shù)據(jù)級并行(SIMD)奠定了基礎(chǔ)。

(2)模板計算與優(yōu)化。

所謂模板計算，指的是多次迭代，且在一次迭代內(nèi)按網(wǎng)格點的順序, 依次對所有網(wǎng)格點進(jìn)行更新操作, 更新時會用到該網(wǎng)格點的相鄰點的信息。模板計算可按維度和更新所用鄰居點的數(shù)目來分類，圖1給出了2D-5點、3D-7點模板示例，結(jié)合前面的敘述可知，3D GVF場計算中對每個像素點梯度向量的更新屬于典型的3D-7點模板計算。

Figure 1 Samples of stencil computations圖1 模板計算示例

模板計算的兩個顯著特點是：①不連續(xù)的內(nèi)存訪問模式，容易造成緩存缺失。文獻(xiàn)[9]對模板計算的緩存缺失因素作了詳細(xì)的分析，概括來說，當(dāng)數(shù)組大于緩存容量時, 本次更新的數(shù)據(jù)在下次更新前已經(jīng)被寫回內(nèi)存；而且當(dāng)數(shù)據(jù)量大時，多次迭代會導(dǎo)致數(shù)據(jù)緩存的容量缺失(CapacityMiss)。②計算/訪存比低，緩存中數(shù)據(jù)重用率低，對訪存帶寬要求高。

對模板計算優(yōu)化的關(guān)鍵在于充分開發(fā)計算和數(shù)據(jù)的時/空局部性。諸多優(yōu)化方法中，緩存分塊(CacheBlocking)是一種典型的優(yōu)化思路。RiveraG等[10]提出的緩存分塊策略如圖2所示，其中I是單元跨度(Unit-stride)維度，或稱變化最快的維度；K是變化最慢的維度。實驗表明：針對I和J兩個維度的分塊打斷了線程內(nèi)存讀取流的持續(xù)性，不利于數(shù)據(jù)并行化的展開，因此探索高效的分塊方式，并從XeonPhi體系結(jié)構(gòu)上找到理論支撐，具有重要的應(yīng)用價值。

Figure 2 Rivera cache blocking圖2 Rivera緩存分塊策略

(3)XeonPhiTMcoprocessors體系結(jié)構(gòu)。

XeonPhi擁有數(shù)十個核，每個核包含一個支持512位SIMD的向量處理單元(VPU)，可同時處理8路雙精度或16路單精度浮點數(shù)據(jù)。Intel提供Intel向量化庫、IntelIntrinsic函數(shù)等方式使用該VPU單元。每個核擁有32KB的一級數(shù)據(jù)緩存和32KB的一級指令緩存。此外，每個核還可以使用512KB的L2級Cache。不同核的L2Cache通過雙向內(nèi)存控制器相連。

3 GVF并行算法設(shè)計

3.1 線程級并行和數(shù)據(jù)級并行

圖像空間中各像素點計算的天然并行性，使得XeonPhi眾核集成架構(gòu)的強大計算能力能得到有效發(fā)揮。本節(jié)討論XeonPhi對圖像處理算法的通用加速方法，主要由兩個層次組成：(1)線程級并行；(2)數(shù)據(jù)級并行。此外，還對XeonPhi圖像處理中的存取模式進(jìn)行了有效的探索，可視為XeonPhi對圖像操作的通用預(yù)處理。

首先是圖像邊緣點的判斷和處理，由第2節(jié)的數(shù)值解法可知，計算邊緣點的laplacian時，其鄰居點有的并不存在，為防止內(nèi)存越界訪問，需要在編碼中引入額外指令判斷當(dāng)前處理的點是否為邊緣點，勢必會影響整體計算效率。為解決該問題，本文對圖像進(jìn)行邊界擴充和鏡像填充，向各個方向均擴充一個像素。為方便觀察，圖3以二維圖像的左上角作示例，3D圖像的處理與此一致。

Figure 3 Expanding 2D image with 1 pixel in all directions， the arrows indicates values the new boundary pixels use.圖3 二維圖像邊界擴充與鏡像填充示例， 箭頭指示數(shù)據(jù)拷貝的方向

為使得XeonPhi的512位VPU可以高效地訪問數(shù)據(jù)，提高計算吞吐率，采取如下措施：(1)動態(tài)分配內(nèi)存時，使用_mm_malloc函數(shù)，確保所分配內(nèi)存邊界對齊；(2)對單元跨度維度，即變化最快的維度，在圖像邊界擴充的基礎(chǔ)上再次擴展，確保其長度為64字節(jié)的整數(shù)倍。

根據(jù)圖像數(shù)據(jù)在內(nèi)存中的組織方式，在變化慢的兩個維度使用線程級并行，使用OPENMP編譯指導(dǎo)語句將任務(wù)劃分給多個核；需要說明的是：(1)對XeonPhi而言，線程綁定方式對于計算效率有較大的影響，使用“KMP_AFFINITY=balanced”模式確保所有線程平衡地劃分到Xeon Phi協(xié)處理器的核上；(2)使用OPENMP提供的collapse指導(dǎo)語句將兩個維度的循環(huán)折疊到一個大的循環(huán)中，能有效減少OPENMP任務(wù)調(diào)度的開銷。

低維度并行方面，使用編譯指導(dǎo)語句確保Xeon Phi 512位的SIMD單元得到有效的利用。

3.2 緩存分塊策略

本節(jié)將模板優(yōu)化相關(guān)方法引入3D GVF場計算加速，在綜合考慮3D GVF模板計算特點和Xeon Phi體系結(jié)構(gòu)的基礎(chǔ)上，提出一種Xeon Phi平臺相關(guān)的分塊策略。

首先，為充分發(fā)揮Xeon Phi中SIMD單元的效能，對單元跨度維度，即變化最快的維度不分塊，如此能保證每個線程的內(nèi)存讀取流持久而連續(xù)；其次，參考文獻(xiàn)[11]，對N×N×N大小的區(qū)域，推薦的形狀為(N-2)×s×(s×L/2)，其中，L為一個Cache行的長度，s為分塊大小。對Xeon Phi而言，每個核的L2緩存為512 KB，L=64B，為確保內(nèi)存讀取各鄰居偏移量項時的局部性，也為確保 Xeon Phi 架構(gòu)下使用 512 KB L2 高速緩存時的局部性，避免Cache容量缺失，近似有：

N×s×(s×64/2)×Tp×Nm<512KB

(3)

3.3 整體算法(偽代碼)

3D GVF算法加速偽代碼如下所示：

1. for (t=0;t

2. #pragma omp parallel for collapse(2)

3. for(jj=1;jj

4. for(kk=1;kk

5. for(k=kk;k

6. for(j=jj;j

7. #pragma simd

9. …// 3D GVF stencil computation

10. }

11. }

12. }

13. }

14. }

15. }

其中,timesteps指迭代步數(shù);nx、ny、nz分別代表圖像的三個維度大小;s為第二維分塊大小，32s為第三維分塊大小。

1.3 利用果蠅的蛹收集處女蠅 接種純種親本果蠅，當(dāng)觀察瓶壁上出現(xiàn)較多的黑褐色的蛹時，用干凈解剖針輕輕地把黑褐色的蛹取出，單獨放到10mL的塑料離心管里，置于溫度25℃、濕度60%的培養(yǎng)箱里進(jìn)行培養(yǎng)。每天觀察，待果蠅羽化出來后進(jìn)行麻醉，鑒別雌雄，收集處女蠅。

4 實驗結(jié)果分析

4.1 實驗平臺

本文以Native模式使用XeonPhi，此時其可視為一個獨立的處理器，完成所有運算，不受CPU控制。此模式下運行的代碼除了512位SIMD指令外，不能含有其他SIMD指令，而且編譯時需要加入編譯選項-mmic。作為對比，本文亦在8核2.60GHz的IntelXeonE5-2670CPU上對3DGVF算法進(jìn)行了性能測試(采用AVX進(jìn)行向量化)。XeonPhi與E5-2670CPU的具體配置如表1所示。采用的編譯器為Inteliccversion13.0，OPENMP版本為3.1。

Table 1 System configuration表1 系統(tǒng)配置

為了克服隨機因素的影響，本文所有測試采用執(zhí)行5遍求平均值的方式進(jìn)行，每遍執(zhí)行時迭代100次，針對不同的問題規(guī)模，采用關(guān)鍵函數(shù)墻內(nèi)時間作為性能度量。

4.2 多線程和SIMD對GVF場計算速度的改進(jìn)

本小節(jié)評測在Xeon Phi上使用多線程和向量化取得的性能提升。圖4是不同線程下，有無向量化時的3D GVF場計算性能對比，其中，橫坐標(biāo)是所用線程數(shù)，縱坐標(biāo)是墻內(nèi)計算時間，單位為秒。

Figure 4 Computation performance under different setups圖4 不同配置下的系統(tǒng)計算性能

可以看出，在眾核集成架構(gòu)中，多線程能顯著提升計算性能，在線程數(shù)較低時表現(xiàn)尤其明顯。但是，當(dāng)線程數(shù)達(dá)到一定規(guī)模時(滿足每核兩個線程時)，線程數(shù)目的提升對計算性能的影響不大。對向量化SIMD而言，當(dāng)線程數(shù)較少時，使用向量化能以超過兩倍的加速比提升計算性能。但是，當(dāng)線程數(shù)較多時，SIMD帶來的性能改進(jìn)已幾乎可忽略不計(<1%)，本次測試中，在線程數(shù)較多時，甚至出現(xiàn)了SIMD拖低性能的現(xiàn)象。綜合而言，使用Xeon Phi進(jìn)行圖像相關(guān)算法處理，尤其是像素級stencil遍歷運算時，滿足每核兩個線程是較優(yōu)的選擇。

4.3 不同分塊

本小節(jié)評測提出的分塊策略能否達(dá)到較優(yōu)的效果。與之前諸多采用遍歷不同分塊求最優(yōu)的工作不同，本文給出了分塊大小的經(jīng)驗指導(dǎo)公式，如公式(3)所示。為評測該公式的性能，設(shè)計如下測試：問題規(guī)模選取5123，將第二維分塊的大小即第3.2節(jié)中的s從1到8遍歷，第三維分塊的大小選取為n×s，n的取值為1、2、4、8、16、32。本文經(jīng)驗公式給出的分塊大小為(512,2,64)。

將Xeon Phi的線程數(shù)設(shè)置為每核兩個線程，總計114個線程。不采用分塊策略時，Xeon Phi計算51233D GVF場所需時間平均值為8.471秒，以此為基準(zhǔn)，不同分塊大小相對此基準(zhǔn)的加速比如圖5所示，兩個坐標(biāo)軸分別對應(yīng)n(橫坐標(biāo))和s(縱坐標(biāo))的大小。

Figure 5 Speedup ratios under different cache blocking sizes圖5 不同分塊大小相比基準(zhǔn)情況的加速比

可以看出，本文所提出的分塊策略相比基準(zhǔn)情況有1.25倍的加速比，在所測試的組合中是最優(yōu)的，從實驗角度驗證了3.2節(jié)中對3D GVF計算和Xeon Phi L2緩存的分析。觀察圖5還可發(fā)現(xiàn),加速比相對較高的組合集中在圖5所示矩陣的反對角線方向。

4.4 Xeon Phi與Xeon CPU優(yōu)化后的計算性能對比

為直觀體現(xiàn)本文在Xeon Phi上對3D GVF場計算的加速效果，本小節(jié)將Xeon CPU和Xeon Phi的計算性能進(jìn)行了比對，兩者所用優(yōu)化措施均包括多線程和向量化，在Xeon Phi端使用本文提出的經(jīng)驗分塊公式，在Xeon CPU端以遍歷尋優(yōu)的方式取最佳性能。兩者的對比結(jié)果如表2所示，可見本文的優(yōu)化方法使得Xeon Phi的絕對性能達(dá)到了Xeon CPU的2.77倍。

Table 2 Performance comparisonsbetween Xeon CPU and Xeon Phi表2 Xeon CPU和Xeon Phi的性能比較

5 結(jié)束語

本文首次針對Xeon Phi平臺進(jìn)行了GVF場計算加速研究。GVF場的計算體現(xiàn)出圖像處理中各像素的天然并行性，在合理安排圖像數(shù)據(jù)內(nèi)存存取模式的基礎(chǔ)上，結(jié)合通用的線程級并行和數(shù)據(jù)級并行對GVF場進(jìn)行加速優(yōu)化；同時，GVF場計算是一種典型的模板計算，結(jié)合Xeon Phi二級Cache結(jié)構(gòu)特點，提出了高效的分塊大小經(jīng)驗公式，避免了費時費力的分塊尋優(yōu)。實驗結(jié)果表明，本文所提方法在Xeon Phi平臺計算3D GVF場取得了很好的加速比。

[1] Kass M, Withkin A, Terzopoulos D. Snakes:Active contour models[J]. International Journal of Computer Vision, 1988,1(4):321-331.

[2] Zhong Y, Jain A K, Dubuisson-Jolly M P. Object tracking using deformable templates[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2000, 22(5):544-549.

[3] Caselles V, Morel J M, Sbert C. An axiomatic approach to image interpolation[J]. IEEE Transactions on Image Processing, 1998, 7(3):376-386.

[4] Xu C, Prince J L. Snakes, shapes, and gradient vector flow[J]. IEEE Transactions on Image Processing, 1998, 7(3):359-369.

[5] Jifeng N,Chengke W,Shigang L,et al.NGVF:An improved external force field for active contour model[J]. Pattern Recognition Letters, 2007, 28(1):58-63.

[6] Wang Y Q, Jia Y D. A novel approach for segmentation of cardiac magnetic resonance images[J]. Chinese Journal of Computers, 2007, 30(1):129-136.(in Chinese)

[7] He Z, Kuester F. GPU-based active contour segmentation using gradient vector flow[M]∥Advances in Visual Computing, Berlin:Springer, 2006:191-201.

[8] Smistad E, Elster A C, Lindseth F. Real-time gradient vector flow on GPUs using OpenCL[J]. Journal of Real-Time Image Processing, 2012,DOI 10.1007/S11554-012-0257-6.

[9] Leopold C. Cache miss analysis of 2D stencil codes with tiled time loop[J]. International Journal of Foundations of Computer Science, 2003, 14(1):39-58.

[10] Rivera G, Tseng C W. Tiling optimizations for 3D scientific computations[C]∥Proc of ACM/IEEE 2000 Conference on Supercomputing, 2000:32.

[11] Leopold C. Tight bounds on capacity misses for 3D stencil codes[C]∥Proc of the International Conference on Computational Science-Part I, 2002:843-852.

附中文參考文獻(xiàn)：

[6] 王元全, 賈云得. 一種新的心臟核磁共振圖像分割方法[J]. 計算機學(xué)報, 2007, 30(1):129-136.

QIJin,born in 1988,MS candidate,his research interest includes system software.

李寬(1984-),男，山東寧陽人，博士，助理研究員，研究方向為并行計算和圖像處理。E-mail:likuan@nudt.edu.cn

LIKuan,born in 1984,PhD,assistant researcher,his research interests include parallel computing, and image processing.

楊燦群(1968-),男，湖南桃江人，博士，研究員，研究方向為系統(tǒng)軟件。E-mail:canqun@nudt.edu.cn

YANGCan-qun,born in 1968,PhD,research fellow,his research interest includes system software.

杜云飛(1980-),男，安徽阜南人，博士，助理研究員，研究方向為并行計算、編譯技術(shù)和程序性能優(yōu)化。E-mail:forest80@163.com

DUYun-fei,born in 1980,PhD,assistant researcher,his research interests include parallel computing,compiler technology, and program performance optimization.

Accelerating3DGVFfieldcomputationonXeonPhiusingstenciloptimization

QI Jin1,LI Kuan2,YANG Can-qun1,DU Yun-fei2

(1.National Laboratory of Parallel and Distributed Processing,National University of Defense Technology,Changsha 410073;(2.College of Computer Science,National University of Defense Technology,Changsha 410073,China)

3D Gradient Vector Flow (GVF) field has wide applications in many image processing algorithms. The computation of GVF field typically needs several iterations and is rather time consuming. Therefore, it is important and meaningful to improve the computation speed of 3D GVF field. The data level parallelism and thread level parallelism are introduced to accelerate the GVF field computation procedure on Intel Xeon Phi many core integrated platform for the first time. Meanwhile, GVF field computation is a kind of stencil computation, whose computation-memory access ratio is low. A novel cache blocking strategy is proposed to fully utilize the L2 cache of Xeon Phi architecture，and to improve the computation speed of GVF field. The experimental results show that the proposed optimizations could effectively improve the speed of GVF filed computation. Especially, for a 51233D image, compared with the performance obtained by a 2.6G Hz 8 core 16threads Intel Xeon E5-2670 CPU, the speedup achieved on Xeon Phi is 2.77X.

3D GVF field;Xeon Phi;stencil optimization;cache blocking

1007-130X(2014)08-1435-06

2013-08-12;

：2013-11-11

國家863計劃資助項目(2012AA010903)；國家自然科學(xué)基金資助項目(61170049，61303189)

TP393

：A

10.3969/j.issn.1007-130X.2014.08.003

齊金(1988-),男，湖南株洲人，碩士生，研究方向為系統(tǒng)軟件。E-mail:qijin2012@yeah.net

通信地址：410073 湖南省長沙市國防科學(xué)技術(shù)大學(xué)并行與分布處理重點實驗室

Address:National Laboratory of Parallel and Distributed Processing,National University of Defense Technology,Changsha 410073,Hunan,P.R.China