基于ＧＰＵ的３Ｄ距離變換
2008-12-31 00:00:00田緒紅司徒志遠陳茂資韓國強
計算機應用研究 2008年9期
摘 要：距離變換在圖像處理中有著非常廣泛的應用。由于3D圖像數據的復雜性，傳統基于CPU的3D距離變換效率較低。為此，研究了將3D圖像數據有效地組織到紋理中存儲的方法，設計并實現了基于GPU的3D距離變換并行算法。實驗結果表明，相對基于CPU的算法，該方法具有非常高的加速比。
關鍵詞：三維距離變換； 圖形處理器通用計算； 并行算法
中圖分類號：TP391.4 文獻標志碼：A
文章編號：10013695(2008)09284703
3D distance transform on GPU
TIAN Xuhong1，2， SITU Zhiyuan1， CHEN Maozi1， HAN Guoqiang
(1. College of Informatics， South China Agricultural University， Guangzhou 510642， China； 2. College of Computer Science Engineering， South China University of Technology， Guangzhou 510642， China)
Abstract:Distance transform had gained numerous applications in image processing. Traditional CPUbased 3D distance transform was inefficient because of the complexity of 3D images. This paper proposed an efficient method to store 3D images in the texture of graphics hardware. A parallel 3D distance transform algorithm based on GPU was implemented. Experiments show that the proposed GPUbased method gain very high speedup ratio compared with CPUbased algorithm.
Key words：3D distance transform; GPGPU; parallel algorithm
0 引言
距離變換是數字圖像處理中的一個經典問題，在骨架抽取、形狀匹配、目標重建、圖像分析與模式識別等許多方面都有非常廣泛的應用。傳統基于CPU的3D距離變換算法效率較低，如一幅大小為N×N×N的3D圖像，常用的算法時間復雜度為O(N4/2)。
近些年來圖形處理器(GPU)的硬件技術飛速發展，平均每年就有新的一代GPU推出市場。與之相應的是其可編程能力不斷提高，出現了如C for graphics (CG)、high level shading language(HLSL)、OpenGL shading language(GLSL)等類似C的高級語言。GPU的應用領域不斷擴大，學者們提出GPU通用計算(GPGPU)的概念，并成功地將GPU應用于非圖像繪制方面的計算。這些計算相當廣泛，如碰撞檢測[1]、數值計算[2]、流體模擬[3]、圖像處理[4]、偏微分方程計算[2，5，6]、3D網格模型細分[7，8]等。由于GPU具有高效的并行處理能力，使其在許多通用計算領域有很好的性能，達到一些大型計算機的水平，并且GPU價格低廉，有利于普及推廣。
距離變換方面，Strzodka等人[9]提出了利用GPU進行2D廣義距離變換；Sigg等人[10]在GPU上實現了2D帶符號的距離變換；Sud等人[11]提出了基于GPU的3D網格模型的距離變換，而基于GPU的3D圖像距離變換尚未有報道。通過對3D圖像距離變換算法及GPU通用計算原理的分析，本文提出一種利用GPU的并行處理能力進行3D圖像距離變換的方法。
1 GPU通用計算
GPU是專門為圖形圖像處理而設計的硬件芯片，其設計初衷是為了減輕CPU的負擔，將大部分的圖形圖像運算如光照處理、頂點運算、紋理映射等功能轉移到一塊專門的芯片上，這就是圖形處理器(GPU)的由來。從系統架構上看，GPU 是針對向量計算進行了優化的高度并行的數據流處理器。其中包括兩種流處理單元：頂點著色器(vertex shader)，是多指令多數據流的處理單元(MIMD)；像素著色器 (pixel shader)，是單指令多數據流的處理單元(SIMD)[12，13]。GPU的系統架構決定了其在數據流處理上可以獲取較高的效率。但GPU在系統架構設計上對整數運算、邏輯運算支持的貧乏，以及其不同于傳統CPU的數據存取方式，使其在非圖形圖像計算方面也存在一些局限。GPU在通用計算方面不是萬能的，它只能完成一些特定的計算。一般說來GPU通用計算有如下特點[13，14]： a)算法程序必須能被分解成為不相關的程序片段(section)，這些程序片段作為GPU的計算內核，每個流處理器執行相同的程序片段來并行求解最終的計算結果；
b)由于GPU的特殊架構，使之更適合運算程序的循環結構部分而不適合運算程序的分支結構部分；
c)每個程序片段的控制結構不能復雜，語句不能太多，實際應用中可以將復雜算法程序分解成幾個規模較小的代碼模塊；
d)由于GPU是針對圖形處理進行加速的，其中最主要的計算元素，或者說數據結構就是紋理。
2 基于GPU的3D距離變換算法
2.1 紋理的組織
如前文所述，GPU通用計算的數據對象主要是紋理。因此，利用GPU的并行計算能力實現3D距離變換算法，首先要解決的一個關鍵問題是底層數據的表示。具體來說就是如何將3D圖像數據有效地組織到紋理中。雖然GPU本身提供3D紋理的操作，但是其對3D紋理的支持尚不成熟，效率也不高，而相對來說由于對2D紋理的渲染本身就是GPU的一個設計目標，其對2D紋理提供了高效的處理能力。
對于一幅N×N×N的3D二值圖像，可以將其分解為N個N×N的2D二值圖像。3D距離變換算法中，對每個前景點尋找其26鄰域點只需要相鄰的3幅2D圖像。從理論上說，雖然可以將3D體數據存放到內存的3D數組中，每次GPU處理時再將需要的數據寫入紋理，最后將計算得到的紋理寫回內存，但是這樣將導致內存和顯存頻繁的數據交換，嚴重影響算法的執行效率。本文使用一個大小為N的2D紋理數組將上面的N個N×N的2D二值圖像全部讀入顯存。一個標準的紋理單元有四個通道(ARGB)，占用4×16＝64 bit的空間。考慮當今主流的顯存容量為256 MB，可以存放256 MB/64 bit=4 194 304個紋理單元，約只能存放322×322×322大小的3D數據。由于距離變換是對二值圖像進行操作，得到的結果是一個灰度圖像，只需一個通道就可以了。可以考慮將四個數據元壓縮進一個紋理單元中，即一個紋理單元的四個通道分別存放一個數據元，這樣256 MB的顯存就能存放下512×512×512大小的3D數據集。當然主流顯卡都能夠設置一部分內存作為顯存的補充，這樣就能放下更大的3D圖像了，但性能上有所損失。
2.2 算法流程
本文采用微軟公司的Direct3D作為應用程序的接口，用其自帶的高級著色語言HLSL編寫了3D距離變換的著色器程序，即流處理器程序片段。該程序片段在Direct3D中又稱為效果。在該效果中實現兩個手法(technique)：a)一個用于實現距離變換的初始化，即將背景距離值設為0，前景點距離值設為無窮大(具體實現時可設置為一個最大可能的距離變換值)；b）距離變換的核心算法，根據第i-1、i、i+1三層紋理確定第i層每一紋素的26鄰域點的距離值，并以紋理渲染的方式更新其距離值。
此處3D距離變換采用Chamfer距離測度。令p、q為3D圖像中相鄰兩點，則它們之間的Chamfer距離定義為
其中：n1=3，n2=4，n3=5。
于是3D Chamfer距離變換定義為
其中：Nm(p)為p的m鄰域，通常取m=26；F表示前景點集合。
整個距離變換算法流程如圖1所示。由圖1可以看出，本算法將運算量最大的部分都放進了GPU中運行，這樣就能充分利用GPU的并行處理能力，以提高算法的效率。
下面以距離變換一次迭代為例，詳細說明其實現過程：
a)依次將第i層及其相鄰的第i-1，i+1層紋理用effect>setTexture()方法設置為著色器程序可以使用的紋理。
b)用effect>setVectorArray()方法將vSampleOffset數組(存放單層紋理中8鄰域的紋理坐標偏移量)傳入著色器。
c)用device>getRenderTarget()方法將一個臨時的紋理tempTexture設置為渲染目標，用于保存計算結果(紋理不能同時讀寫，故需要一個臨時紋理用于寫入)。
d)根據Chamfer距離變換定義，在GPU中更新第i層紋理的每個像素點的距離值，結果保存在臨時紋理中，如圖2所示。這一步驟使用的HLSL代碼如下：
PS_OUTPUT psDT(VS_INOUTPUT In)
{
PS_OUTPUT Output; /* 聲明一個存儲返回值的結構體變量 */
float v1， v2， v3;
vSampleOffset[0].xy += In.vTexCoord.xy;
……
vSampleOffset[8].xy += In.vTexCoord.xy;
/* 從CPU傳入的8鄰域紋理坐標的偏移量 */
v1 = min3(
tex2D(sImage0， vSampleOffset[0].xy).r + (PACE*5)，
tex2D(sImage1， vSampleOffset[0].xy).r + (PACE*4)，
tex2D(sImage2， vSampleOffset[0].xy).r + (PACE*5));
v2 = min3(
tex2D(sImage0， vSampleOffset[1].xy).r + (PACE*4)，
tex2D(sImage1， vSampleOffset[1].xy).r + (PACE*3)，
tex2D(sImage2， vSampleOffset[1].xy).r + (PACE*4));
v3 = min3(
tex2D(sImage0， vSampleOffset[2].xy).r + (PACE*5)，
tex2D(sImage1， vSampleOffset[2].xy).r + (PACE*4)，
tex2D(sImage2， vSampleOffset[2].xy).r + (PACE*5));
……
v1 = min3(v1， v2， v3);
/* 根據當前點的26鄰域點的距離值，更新該點的距離值 */
Output.vColor = float4(v1， 0， 0， 0);
return Output; /* 像素著色器返回值，用于渲染到場景 */
}
e)將渲染結果tempTexture寫回第i層紋理，對第i+1層紋理重復以上步驟。
3 實驗結果
實驗采用的數據如圖3和4所示。其中圖3為一組形狀相同但大小不等的規則立方體圖像，全部由前景點構成；圖4為非規則的實際3D重建根系圖像，其前景點比較稀疏。
實驗平臺為：2.0 GHz的AMD雙核速龍3600+處理器，1 GB 雙通道DDR2667內存，Windows XP SP2操作系統。GPU為NVIDIA公司的GeForce 7600GT，其核心工作頻率為575 MHz，顯存頻率為1 400 MHz，顯存位寬為128 bit，顯存大小為256 MB。它提供了8個頂點處理器及12個像素處理器，每次能并行處理含4位通道的1個頂點或像素。應用程序使用C語言進行編寫，采用Direct3D 9c API以及HLSL來進行GPU方面的程序設計。在Visual Studio 2005的開發環境下進行了程序的編譯及執行。
為了與CPU計算的效率進行對比，對上述數據采用基于CPU的算法進行了實驗。圖5~7給出了在CPU與GPU上進行距離變換的計算時間比較。從圖中可以看出，相對CPU而言，GPU進行距離變換的速度得到明顯提高。
初始化過程GPU所用時間非常少，而在CPU上運行的時間則與圖像大小成正比。
距離變換迭代階段，對于圖3的規則圖像，CPU運行時間符合時間復雜度O(N4/2)；對于圖4非規則圖像，由于有許多背景點，這些點在距離變換過程中是不用處理的，因而CPU運行時間正比于M×I。其中：M為圖像前景點數；I為最外層循環迭代次數，由圖像中前景物體的最大橫截面半徑決定。
GPU進行距離變換，理論上加速比最大不超過像素處理器的個數。本文實驗平臺像素處理器的個數為12，因此，理論上加速比最大為12。但由圖7可以看出，當數據量達到一定規模后，實際加速比遠遠超過12。這是由于GPU 是針對向量計算進行了優化的高度并行的數據流處理器，GPU計算時，紋理是作為整體加載到流處理單元中執行，因而不能簡單地用類似 CPU中的時間復雜度推算GPU的運行時間。
由圖7也可發現，當數據量較小時，GPU加速比明顯降低。其主要原因是數據量小時，流處理器的優勢沒有發揮出來，而紋理的載入與載出有一定的消耗，使整體加速效果下降。但由于GPU初始化時間非常少，即使是前景點相對較少的實際根系圖像，GPU總加速比仍然達到4.6以上。
圖8、9為GPU距離變換的效果圖。其中，圖8為圖3（b)立方體距離變換后的結果切掉1/4展示的內部效果圖。圖9為圖4（b)橫切Y軸的連續16個截面圖。效果圖表明GPU距離變換的結果是正確的。
4 討論
近年來也有不少學者研究CPU上距離變換的改進算法，目前有些改進算法的時間復雜度可提高到O(N3 log2n）。其中：N3是圖像大小；log2n是外層循環的次數。對于256×256×256的規則圖像，這種改進算法的加速比為（256/2）/log2 256=16，遠遠小于GPU上算法加速比。當然隨著圖像大小增加，這種CPU上的改進算法加速比會提高，但也遠小于GPU上的算法。而且更重要的是對于實際3D圖像，外層循環只由前景物體的最大橫截面半徑決定，通常遠遠小于前景物體所處的圖像空間的大小，因而，CPU上改進的距離變換算法效率提高相對較小。
參考文獻：
［1］GOVINDARAJU N K， REDON S， LIN M， et al. Interactivecollisiondetectionbetweencomplexmodelsin large environmentsusinggraphics hardware[C]//Proc of Eurographics/SIGGRAPH Workshop on Graphics Hardware. 2003:2532.
[2]KRUGER J， WESTERMANN R. Linear algebra operators for GPU implementation of numerical algorithms[J].ACM Trans on Graphics，2003，22(3):908916.
[3]HARRIS M J. Fast fluid dynamics simulation on the GPU[M]//FERNANDO R.GPU Gems Chapter 38. New York: AddisonWesley ，2004.
[4]MORELAND K， ANGEL A. The FFT on a GPU[C]//Proc of the Graphics Hardware. 2003.
[5]BOLZ J， FARMER I， GRINSPUN E， et al. Sparse matrix solvers on the GPU:conjugate gradients and multigrid[J].ACM Trans on Graphics，2003，22(3):917924.
[6]GOODNIGHT N，WOOLLEY C， LUEBKE D，et al. A multigrid solver for boundary value problems using programmable graphics hardware[C]//Proc of the Graphics Hardware.2003:102111.
[7]SHINE L J， JONES I， PETERS J. A realtime subdivision kernel[J].ACM Trans on Graphics，2005，22(3):10101015.
[8]COHENOR D， SLAVIK P. GPUbased composite subdivision[J].Eurographics，2007，26(3):276285.
[9]STRZODKA R， TELEA A. Generalized distance transform and skeletons in graphics hardware[C]//Proc of EG/IEEE TCVG Symposium on Visualization. 2004:221230.
[10]SIGG C， PEIKERT R， GROSS M. Signed distance transform using graphics hardware[C]//Proc of IEEE Visualization.2004:8390.
[11]SUD A， OTADUY M A， MANOCHA D. DiFi:fast 3D distance field computation using graphics hardware[J].Computer Graphics Forum，2004，23(3):557566.
[12]吳恩華. 圖形處理器用于通用計算的技術、現狀及其挑戰[J]. 軟件學報， 2004，15(10):14931504.
[13]OWENS J D， LUEBKE D， GOVINDARAJU N， et al. A survey of generalpurpose computation on graphics hardware[J].Computer Graphics Forum，2007， 26(1):80113.
[14]黃鑫， 李勝， 程惠閣， 等. 基于GPU 的B樣條曲面加速計算[J]. 系統仿真學報，2006，18(增刊1):14.
