基于OpenCL的DRR算法優化研究

田林琳，李 瑩

(沈陽工學院 信息與控制學院，遼寧 撫順 113122)

1 概 述

放射治療計劃系統(TPS)是利用計算機軟件技術、劑量計算方法、控制手段等，解決治療過程中的精確定位、劑量計算的快速精確、保證目標劑量的準確實施以及治療計劃的優化等問題[1]。數字重建放射影像(DRR)是TPS系統中的關鍵功能，是通過模擬X光線的衰減過程，從CT斷層數據中生成類似X光片的成像效果，這對于放療過程中的質量控制具有非常重要的意義[2]。DRR可在任意三維空間角度下實現對病患的“透視”模擬，能夠隨意觀察放射靶區、特定組織或器官，或靶區與周圍器官間的空間關系，因為沒有治療床的限制，可利用DRR得到模擬定位機難以拍攝的圖像。同時省略了使用模擬放射攝影圖像的成像步驟，既節省了臨床費用，也減少了對病患的輻射[3]。

在TPS的臨床應用中，DRR具有重要作用，為了達到臨床應用的目標，要求DRR能夠準確、快速地反映放射治療中的真實情況，并能夠滿足動態性、交互性的需要[4]。這也是三維圖像可視化領域中，人們一直試圖解決的兩個問題：一是實時地進行圖像繪制，在DRR的實際使用中，要求可以快速、多角度地生成圖像；二是真實地顯示結果，這要求計算結果能夠真實、客觀地反映其解剖結構在三維空間上的細節與特征[5]。

典型的DRR算法是以光線跟蹤算法為基礎，利用虛擬攝影機視圖追蹤光線至投影圖像中的每一個像素，每一條光線所經過的體元數據進行累加，累加值與像素亮度值成正比，使用這種方法實現重建出的DRR圖像顯示精度和圖像質量都非常好，滿足了顯示結果真實性的要求[6-7]。但由于此算法需要跟蹤每條從虛擬光源發出的光線，涉及到大量的光線與CT三維體數據進行求交測試的運算，運算量非常大。傳統的基于CPU技術的實現方式通常不能滿足交互式DRR算法的需求，隨著GPU(graphic processing unit)在并行計算能力、存儲容量和可編程能力等方面的發展，使DRR算法的性能提升成為可能[8]。

英偉達(NVIDIA)公司于2007年提出的CUDA架構(compute unified device architecture)使得開發者能夠以GPU的強大并行計算能力為基礎，建立一種更快更高效的數據計算解決方案。但CUDA架構的一個最主要問題是通用性，人們只能使用NVIDIA系列產品對GPU進行開發。而作為一個軟件系統，TPS是運行在醫院提供的計算機上，各醫院的硬件環境千差萬別，所以使用CUDA加速的局限性比較高。而基于異構平臺的OpenCL通用計算規范，允許開發者以相同的代碼在任何支持該規范的平臺上運行，包括NVIDIA和AMD顯卡，甚至包括支持OpenCL的CPU，這大大降低了程序對硬件平臺的依賴性，具有較強的理論意義和實踐價值。

使用OpenCL對TPS系統中的DRR算法進行并行加速，主要出于三個目的：

(1)通過對DRR算法的加速，滿足算法交互性和實時性的需要，這不僅能減少患者的等待時間，也提高了TPS的使用效率。

(2)OpenCL作為通用計算標準，與CUDA架構相比，其性能缺乏評判，通過使用OpenCL技術在NVIDIA顯卡上實現對DRR算法的加速，軟件人員可以更好地審視其性能，以便于在今后的開發中進行平臺選擇。

(3)通過將數據存儲在不同的硬件存儲器上的性能對比，開發人員可以深入了解不同存儲器在并行優化中的影響，便于在使用OpenCL技術時進行數據的合理布局。

2 OpenCL概述

OpenCL是一個可以在異構平臺上進行并行編程的開放框架標準，包括一組底層程序接口和一個高效、快速、可移植的抽象層，為開發者提供了一個有效的并行開發環境、一組與平臺無關的工具和豐富的中間層。

OpenCL框架運行時使用主機端程序對所有支持OpenCL的設備進行統一管理。對于每個支持OpenCL規范的計算設備，其在內部又進一步劃分為多個計算單元(compute unit)，每個計算單元又劃分為多個處理單元(processing element)，每個處理單元以SIMD或SPMD的方式運行[9]。

為了支持代碼的可移植性，OpenCL定義了一個抽象的內存模型，開發者可以根據該模型編寫代碼，硬件廠商可以將該模型映射到他們自己實際的硬件上[10-11]。OpenCL定義的內存空間如圖1所示，這些內存空間與OpenCL程序密切相關。

圖1 OpenCL內存模型

OpenCL內存模型按照訪問權限劃分，依次分為全局內存和常量內存(Kernel范圍)、本地內存(Workgroup范圍)、私有內存(Work-item范圍)。其中紋理內存和常量內存的數據位于全局內存，但其擁有高速緩存[12]。

3 DRR算法的優化過程

3.1 DRR算法原理

DRR圖像是通過模擬X光線衰減過程生成的，X光線穿過的組織不同，被吸收的劑量也不同，所以采用X光線衰減模型進行模擬，有如下公式：

(1)

其中，I0為X射線的初始強度；μi為組織i的線性衰減系數；L為X射線的長度。

式1的離散形式如式2所示：

I=I0×e-∑μili

(2)

其中，I0和μi同式1；li表示X射線穿過組織i的距離。式2為生成DRR圖像的關鍵公式[13]。

3.2 編寫串行程序

編寫一個能夠正確運行、沒有存儲器資源泄露的CPU串行程序，作為并行優化的正確性和性能測試的基準。這樣做的原因是：(1)由于在目前的OpenCL版本上，使用C++只能實現一些OpenCL編程中主機部分的函數，包括OpenCL設備初始化、CPU和GPU之間的數據通信、注銷OpenCL計算環境等，而真正在GPU上執行的kernel函數和device函數都是以C語言函數的形式編寫的，所以首先把TPS系統中DRR算法部分由C++程序轉換為C程序。(2)OpenCL device上運行的代碼不易調試，所以在串行程序中去除不易并行的因素，使得程序中的循環或者迭代之間沒有或者較少的數據依賴，這會大大減少并行程序出錯的概率，因此接下來在第1步串行程序的基礎上改寫以符合OpenCL的并行操作，例如函數中的靜態變量和全局變量在程序運行過程中，會被多個線程同時訪問，但其值一直保持不變，會使并行程序中線程之間存在依賴性，不利于程序的并行化，所以根據算法的實現過程改造為局部變量。還有將程序中代表CT體數據的數據指針由二級指針修改為一級指針，這樣CT體數據成為一個連續的內存數據塊，方便了數據從CPU端向設備端的復制和在綁定紋理時對體數據內存排列的要求。

3.3 確定并行化方案

支持OpenCL計算的GPU通常具有數個多處理器，每個多處理器都具有SIMD架構，支持在同一時鐘周期內處理不同的數據，這使得OpenCL主要針對基于數據的并行計算方式[14]。光線追蹤算法中每一條光線的計算過程都是相互獨立的，不同光線之間沒有依賴關系，因此可以方便地進行并行化，只需將每一條光線設計成一個獨立的計算線程就能實現基于光線追蹤的DRR算法并行化。這是利用OpenCL實現基于光線追蹤的DRR算法的可能性。圖2是DRR算法的并行化方案。

圖2 DRR算法的并行化方案

使用OpenCL進行并行程序設計，其中一個主要特點是針對不同的存儲器類型進行數據存儲。GPU中有全局內存、常量內存、紋理內存和共享內存。DRR算法的輸入數據CT體數據少則80～150張，多則200～400張甚至更多。以150張CT圖像為例，其需要的內存為150*0.5 MB=75 MB，所以只有全局內存和紋理內存可供選擇。紋理內存緩存在芯片上，因此它能夠減少對內存的請求并提供更高效的內存帶寬[15]。紋理緩存是專門為那些在內存訪問模式中存在大量空間局部性(spatial locality)的圖形應用程序而設計的[16]。文中對CT體數據存儲在全局內存和紋理內存的性能進行了對比。

常量內存中的數據位于全局內存，但擁有局部緩存加速，常量內存的空間較小(只有64 K)，因此將虛擬光源的點坐標、CT圖像的層厚、圖像大小、像素大小等在計算過程中恒定不變的少量參數放在常量內存中。

對于不同CT值數據的衰減表，大小是4 096，類型是單精度浮點，可以選擇將數據存儲到全局內存、紋理內存或者常量內存，文中也對三種情況下的程序性能進行了對比。

3.4 并行化實現步驟

OpenCL實現DRR算法的步驟主要有：

(1)調用clGetPlatformIDs獲得可用的平臺，設置使用OpenCL設備；

(2)調用clCreateContext在指定平臺上創建上下文；

(3)調用clCreateCommandQueue函數，在已經創建好的上下文上創建命令隊列；

(4)在GPU上分配三維image對象，并將CT體數據和CT衰減表從CPU讀取至設備的全局內存或者紋理內存中；

(5)在GPU上為算法的其他參數分配設備內存或者常量設備內存，如虛擬光源的點坐標、CT圖像的層厚、圖像大小、像素大小等，并將CPU端的數據復制到對應的設備內存中；

(6)創建一個網格，并且設置好workgroup和workitem的大小；

(7)創建內核對象，使用clCreateKernel創建一個kernel對象，調用clEnqueueNDRangeKernel啟動kernel函數；

(8)kernel函數執行，先判斷當前的workitem要處理的像素是否在DRR圖像范圍內，如果不是，則不進行計算直接返回；否則執行DRR算法內核函數，在內核函數的執行過程中更新DRR圖像關聯的全局設備內存；

(9)kernel函數執行完畢，把更新后全局內存數據復制回CPU內存；

(10)在屏幕上描畫CPU內存上的結果；

(11)釋放OpenCL設備上的空間以及CPU上的內存。

3.5 并行參數設置

基于OpenCL的DRR算法中，用一個線程來計算DRR圖像中一個像素的值，如DRR圖像的分辨率為m*n，則在一個kernel中就包括使用m*n個線程，每個workgroup中的workitem設置為32*16，則一個kernel中的workgroup數目為(m/32)*(n/16)。

4 實 驗

4.1 實驗環境

針對相同的算法復雜度，將CT體數據分別存儲在全局內存和紋理內存并分別與串行程序進行了對比，另外對CT衰減表分別存儲于全局內存、紋理內存和常數內存并與串行程序進行了對比，以評估OpenCL將數據存儲在不同存儲器上的性能。

使用的平臺如下：

GPU：NVIDIA GeForce GTX760 2 G顯存；

CPU：Intel i5 4590 4核4線程 3.3 GHz；

內存：DDR3 1 600 MHz 8 G。

4.2 結果分析

文中分別對分辨率為2 048*2 048和1 024*1 024的DDR圖像進行了CPU和GPU的計算結果和性能對比，使用GPU運算的結果和CPU相比，運行結果誤差是1e-2，結果可以接受，如圖3所示。

圖3 CPU和GPU的運算結果對比

性能對比如表1～3所示。

表1 體數據在全局內存時的性能對比

表2 體數據在紋理內存時的性能對比

表3 CT衰減表在不同存儲器上的性能對比

由表1～3中的對比數據可見：經過OpenCL優化的DRR算法速度比CPU串行版本明顯快很多；在OpenCL編程中使用紋理內存可以大大提高程序的加速比；在相同條件下，紋理內存的速度略大于常數內存，并且遠大于全局內存；OpenCL并行程序和串行程序的加速比與基于CUDA的DRR圖像生成算法相近[13]。

5 結束語

通過利用DRR算法的可并行性，使用OpenCL技術對基于光線追蹤的DRR算法進行優化，并在GPU上進行了并行算法調優；通過對OpenCL中的各種存儲器的優化使用，顯著提高了程序的運行效率。目前，OpenCL的使用仍然有很多局限，在諸如復雜的分支和邏輯判斷等方面仍然不如CPU靈活，但隨著GPU的快速發展，相信OpenCL的應用范圍將更加廣闊。

參考文獻:

[1] 閆 鋒.數字重建放射影像生成方法及其應用研究[D].合肥：合肥工業大學，2010.

[2] 戴建榮，胡逸民.圖像引導放療的實現方式[J].中華放射腫瘤學雜志，2006,15(2):132-135.

[3] 吳宜燦，李國麗，陶聲祥，等.精確放射治療系統ARTS的研究與發展[J].中國醫學物理學雜志，2005,22(6):683-690.

[4] FERNANDO R,KILGARD M J.The Cg tutorial:the definitive guide to programmable real-time graphics[M].[s.l.]:Addison-Wesley,2003.

[5] 汪 俊，吳章文，張從華，等.基于射線追蹤的數字重建影像增強技術[J].生物醫學工程學雜志,2005,22(1):125-128.

[6] HUANG Jianbing,CARTER M B.Interactive transparency rendering for large CAD models[J].IEEE Transactions on Visualization and Computer Graphics,2005,11(5):584-595.

[7] LEVORY M.Volume rendering by adaptive refinement[J].Visual Computer,1990,6(1):2-7.

[8] RUSSAKOFF D B,ROHLFING T,MORI K,et al.Fast generation of digitally reconstructed radiographs using attenuation fields with application to 2D-3D image registration[J].IEEE Transactions on Medical Imaging,2005,24(11):1441-1454.

[9] 李 森，李新亮，王 龍，等.基于OpenCL的并行方腔流加速性能分析[J].計算機應用研究,2011,28(4):1401-1403.

[10] 賈海鵬，張云泉，龍國平，等.基于OpenCL的拉普拉斯圖像增強算法優化研究[J].計算機科學,2012,39(5):271-277.

[11] NVIDIA Corporation. NVIDIA OpenCL JumpStart guide[M]．America:NVIDIA,2009．

[12] MUNSHI A．The OpenCL specification[S]．America:Khronos OpenCL Working Group，2009.

[13] 杜曉剛，黨建武，王陽萍.基于CUDA的數字重建影像生成算法[J].計算機科學,2015,42(2):301-305.

[14] 劉 磊.基于GPU的醫學圖像三維重建及可視化技術研究[D].廣州:南方醫科大學,2008.

[15] 劉 鵬，高 軍，雷勛祖,等.一種基于光場的數字重建影像快速生成算法[J]. 南方醫科大學學報,2007,27(10):1537-1539.

[16] NAGY Z,KLEIN R.Depth-peeling for texture-based volume rendering[C]//Proceedings of the 11th Pacific conference on computer graphics and applications.Washington,DC,USA:IEEE Computer Society,2003:429-433.