基于GPU加速和矩陣優化的醫學圖像重建

周　偉　董　鵬　孫業全　李耀武

1(濰坊醫學院設備科　山東 濰坊 261053)

2(濰坊醫學院醫學影像學系　山東 濰坊 261053)

?

基于GPU加速和矩陣優化的醫學圖像重建

周偉1董鵬2孫業全2李耀武2

1(濰坊醫學院設備科山東 濰坊 261053)

2(濰坊醫學院醫學影像學系山東 濰坊 261053)

摘要計算機斷層掃描CT(Computed Tomography)等醫學影像設備產生的投影數據巨大，因此由投影數據重建原始圖像是一個需要大量計算的耗時過程。為了提高CT圖像重建速度，在主流的個人計算機平臺上，將基于GPU(Graphic Processing Unit)并行計算的CUDA(Compute Unified Device Architecture)技術應用于Cimmino同步迭代重建算法，并針對GPU并行架構的特點，使用CSR-ELL聯合存儲格式對數據的存儲進行了優化以進一步提高并行效率。實驗結果顯示，該GPU并行Cimmino重建方法能夠減少72%的圖像重建時間，在保證成像質量的前提下提高了重建速度。

關鍵詞圖像重建CUDA計算機斷層掃描稀疏矩陣

0引言

醫學圖像重建[1]是由物體的檢測數據重建物體圖像的過程，常見的醫學成像設備有CT、MRI、PET等，其中CT作為一種X射線影像設備有著廣泛的應用。CT的成像過程主要分投影和重建兩個階段。為了得到足夠精確的圖像，投影階段需要從不同角度對物體進行投影，并產生重建階段需要的大量復雜的投影數據。在重建階段，通過對投影數據進行數學運算而得到最終的灰度圖像。CT圖像重建方法主要有以FBP[2]為代表的解析法和以Cimmino[3]為代表的迭代法兩類。相比于解析法，以Cimmino作為基礎算法的迭代類方法具有成像質量高、適應性強(適合不完全投影重建情況)等特點，SART、CAV等方法均由Cimmino發展而來[4]。Cimmino算法計算量大、耗時長，在一定程度上制約了該重建方法的應用。隨著計算機技術的發展進步，GPU并行計算成為了一種有效的科學計算加速手段，可以提供高性能的SIMD(Single Instruction Multiple Data)并行計算能力，并已經在醫學圖像處理領域得到了應用[5]。在醫學圖像重建方面，Scherl[6]等人提出了基于解析重建方法的GPU并行算法，取得了優于串行算法的重建速度，但提升有限；Elble[7]等人將GPU并行計算應用到迭代重建法，在專業的Tesla并行計算機上提升了圖像重建速度。

近兩年，由于圖形芯片廠商的重視，GPU的通用計算能力得到了顯著增強，相關的編程開發工具也更加完善[8]。在GPU通用計算領域，英偉達(Nvidia)公司的圖形顯示芯片和CUDA編程工具套件占據了領導地位，任何一臺裝備了英偉達顯示卡的計算機配合CUDA就可進行并行程序設計和運算[9]。該公司最新的基于Kepler架構的產品線當中，作為中低端產品的GK107顯示核心也提供了384個流處理單元，使主流的個人計算機平臺具有了很強的并行計算能力，大大降低了并行計算的應用成本[10]。本文在主流的個人計算機平臺上實現了圖像重建模擬的投影和重建過程，通過應用CUDA技術設計實現了基于GPU并行的Cimmino重建算法，在設計中采用了CSR(Compressed Sparse Row)和ELL聯合的存儲優化方法進一步提高了存儲和運算效率。為了驗證并行Cimmino算法的效率，在裝備了Kepler架構圖形芯片的個人計算機上進行了測試，在不降低成像質量的情況下，并行算法能夠顯著提高重建速度。

1圖像重建原理

1.1圖像重建數據采集

數據采集是CT圖像重建過程的第一個階段，在這個階段X線球管從不同角度對被檢物體發出X線，X線穿過被檢物體后到達探測器形成投影數據。本文采用計算機模擬的方法來產生投影數據，投影數據由以下公式獲得：

Ax=b

(1)

其中，A為系統矩陣，x為被檢測物體，b為投影數據。系統矩陣A中的每一個元素代表一條射線穿過檢測區域某一網格的貢獻值，本文取網格對射線的截距作為貢獻值，如圖1所示。圖中系統矩陣A的某一行元素代表某一條射線穿過檢測區域的各個網格時分別取得的貢獻值，對于正方形檢測區域來說，一條射線最多穿過2×檢測區域寬度-1個網格，導致系統矩陣A中多數元素為零值，因此系統矩陣A為稀疏矩陣。將系統矩陣A與需要重建的圖像x相乘即可得到投影數據b。

圖1　系統矩陣示意圖

1.2Cimmino重建算法

圖像重建是在已知系統矩陣A和投影數據b的情況下，求得被檢測物體灰度圖像x的過程，如式(1)。在眾多同步迭代類圖像重建算法中，Cimmino算法占有重要地位，許多同步迭代類算法均由該算法改進而來。Cimmino算法的計算過程是，首先將一個初始解投向線性方程組(如式(1))的每一行所代表N維空間的超平面，然后求出所有投影點的中心點，再以該中心點為初始點投向所有超平面，以此方式迭代執行，直到完成迭代次數或者滿足某一結束條件為止[11]。以下是Cimmino算法的矩陣表示形式：

xk+1=xk+λATM(b-Axk)

(2)

其中，xk+1為列向量，是經過k+1次迭代后得到的重建圖像，λ為松弛系數，AT為系統矩陣A的轉置矩陣，b為投影數據，對角矩陣M由以下公式求得：

(3)

其中，Ai,j代表系統矩陣A中第i行j列的元素。

2基于CUDA并行計算的Cimmino算法

2.1串行Cimmino算法

串行的Cimmino算法如算法1所示，算法外層為一個循環結構，循環次數為算法的最大迭代次數，算法內部為每一次迭代的計算過程。串行算法中包含大量的矩陣相乘計算，矩陣相乘使用循環結構實現。由于圖像重建數據量大，矩陣維數很高，因此以串行方式執行耗時很長。

算法1串行Cimmino算法

(1) for k ← 0:K

//K為最大迭代次數

(2) M_dp_rxk ← M(b - Axk)

(3) xk+1← xk+ λATM_dp_rxk

(4) end

2.2并行的Cimmino算法

并行的Cimmino算法通過使用CUDA并行編程，利用GPU的多線程計算能力加速矩陣相乘的計算過程。并行算法外層是用來進行迭代計算的循環結構，循環內部包含兩個CUDA內核函數的調用。假設在numAngles個角度對物體進行投影，每個角度使用numRays條射線，重建圖像分辨率為n×n，則并行算法在CUDA調用1中使用numAngles×numRays個GPU線程參與計算，每個線程負責計算M與(b-Axk)相乘的結果列向量M_dp_rxk中的一個元素。在CUDA調用2中，使用n×n個GPU線程計算保存重建圖像數據的xk+1列向量，每個線程負責其中一個維度的求解。

算法2并行Cimmino算法

(1) for k ← 0:K

//K為最大迭代次數

(2) //CUDA調用1，計算M(b - Axk)的結果M_dp_rxk

(3) tid1 ← cuda線程編號

(4) 計算M_dp_rxk[tid1]

//每個 GPU線程負責一個元素的計算

(5) //CUDA調用2，計算xk+1

(6) tid2 ← cuda線程編號

(7) 計算xk+1[tid2]

//每個GPU線程負責一個元素的計算

(8) end

2.3CSR-ELL存儲結構優化

圖像重建過程數據量巨大，涉及大量的矩陣運算，對計算機內存儲器的容量要求很高。對于200×200分辨率圖像的重建，如果投影數據來自179個角度，每個角度179條射線，則系統矩陣將達到40 000行、32 041列，若使用四個字節的單精度浮點數存儲數據，系統矩陣將占用接近5 GB的內存，遠遠超過了主流個人計算機的存儲能力。針對系統矩陣是稀疏矩陣這一特點，本文使用稀疏行壓縮(CSR)格式來存儲系統矩陣。如圖2(b)所示，CSR稀疏矩陣包括Data、Col_Idx和Row_Ptr三個數組。Data數組用于存儲系統矩陣中的非零元素，維度等于系統矩陣中非零元素的個數。Col_Idx數組的維度與Data數組相同，用于存儲Data數組中每一個元素在系統矩陣中的列號。Row_Ptr數組用于存儲系統矩陣中每一行的首個非零元素在Data數組中的索引號。通過使用CSR格式存儲系統矩陣，大大降低了圖像重建對計算機存儲器容量的要求。

圖2　系統矩陣CSR-ELL存儲優化

本文在使用CSR格式存儲系統矩陣的基礎上，結合GPU并行架構的特點，進一步用ELL格式思想優化存儲結構，使優化后的存儲結構能更高效地在GPU上進行并行計算。ELL格式的主要思想是對稀疏矩陣進行填充和轉置以便于矩陣的訪問[12]。如圖2(c)所示，本文首先在存儲順序上將按行優先排列的Data數組轉換為按列排放，實現了多個線程之間的合并訪存。假設有線程t1、t2和t3，三個線程分別處理系統矩陣不同行的數據，在進行矩陣相乘運算時，當t1線程訪問系統矩陣首行的數據1時(如圖2(c))，并行的t2、t3線程需要訪問數據3和5，因為數據1、3和5在經過按列排放優化后的Data數組中相鄰，多次訪存操作可以合并為一次操作[13]，所以在GPU并行時能夠取得更快的存儲訪問速度。其次，為了盡可能使各并行線程具有相同的控制流，本文在Data數組中的相應位置補充0值元素，在不影響計算結果的情況下使稀疏矩陣的各行具有相同的元素數量，這樣各GPU線程可以分配到相同數量的數據，避免了控制流分支，進一步提高了GPU并行的效率。

3實驗方法和實驗結果

3.1實驗方法

為了測試本文提出的基于GPU并行和CSR-ELL存儲優化Cimmino算法的運算效率，以Shepp-Logan頭部模型作為實驗用例對其進行了圖像重建，分析了該算法與非優化并行算法和串行算法在運算效率和存儲效率上的差異。本文實驗使用的硬件設備為一臺主流的便攜式個人計算機，配置了I3-3110M處理器、具有384個流處理單元的GeForce GT740M圖形顯示卡和4 GB內存，程序設計使用CUDA 5.5 Toolkit開發包和C程序語言。在對比實驗中，使用了三種不同分辨率對圖像進行重建，重建需要的投影數據通過來自179個角度以及每個角度179條射線模擬產生，每種測試條件設置下進行5次實驗并取平均值作為最終結果。

3.2實驗結果

為了驗證本文所提方法的可靠性，與串行算法對比了在多種條件設置下的圖像重建結果，結果證明本文實現的基于GPU并行和CSR-ELL存儲優化Cimmino重建方法能夠獲得與串行方法一致的重建圖像。圖3為使用GPU并行算法對Shepp-Logan頭部模型進行重建的圖像。

圖3　Shepp-Logan頭部模型重建，從左至右分別為迭代10、20和150次的重建圖像

為了測試本文提出的基于GPU并行和CSR-ELL存儲優化Cimmino算法的運算效率，對不同分辨率的Shepp-Logan頭部模型進行了圖像重建并進行了分析。從表1的測試結果可以看出，本文提出的CSR-ELL存儲優化算法在運算效率上明顯優于串行和未經存儲優化的GPU并行算法，在實驗中能夠取得接近3.6倍的加速比，節省了72%的重建時間，有效提高了圖像重建速度。在GPU全局內存使用效率方面，為了避免控制流分支而對矩陣進行的補足對齊操作導致CSR-ELL優化算法需要占用更多的存儲器。

表1　實驗結果，表中“時間”列為100次的迭代時間，

4結語

本文將CUDA并行技術與Cimmino圖像重建算法相結合，在主流個人計算機平臺上設計和實現了基于GPU并行的Cimmino算法。在并行算法設計中，使用CSR-ELL聯合存儲格式對數據的存儲進行了優化。實驗結果顯示，本文提出的GPU并行算法在運行效率上明顯優于串行算法，能夠節省72%的重建時間，加速比接近3.6倍，同時圖像的重建質量沒有下降。在今后的研究中，工作將主要集中在如何進一步提高重建速度和優化高分辨率圖像重建的存儲器使用效率方面。

參考文獻

[1] Zeng G L. Medical Image Reconstruction[M]. Heidelberg: Springer, 2010.

[2] Shepp L A, Logan B F. The Fourier reconstruction of a head section[J].Nuclear Science, IEEE Transactions on,1974,21(3):21-43.

[3] Cimmino G. Calcolo approssimato per le soluzioni dei sistemi di equazioni lineari[M].Istituto per le applicazioni del calcolo,1938.

[4] Gilbert P. Iterative methods for the three-dimensional reconstruction of an object from projections[J].Journal of Theoretical Biology,1972,36(1):105-117.

[5] Eklund A, Dufort P, Forsberg D, et al. Medical image processing on the GPU-Past, present and future[J].Medical Image Analysis,2013,17(8):1073-1094.

[6] Scherl H, Keck B, Kowarschik M, et al. Fast GPU-based CT reconstruction using the common unified device architecture (CUDA)[C]//Nuclear Science Symposium Conference Record,2007.NSS’07.IEEE.IEEE,2007,6:4464-4466.

[7] Elble J M, Sahinidis N V, Vouzis P. GPU computing with Kaczmarz’s and other iterative algorithms for linear systems[J].Parallel Computing,2010,36(5):215-231.

[8] Nickolls J, Dally W J. The GPU computing era[J].IEEE Micro,2010,30(2):56-69.

[9] Nickolls J, Buck I, Garland M,et al.Scalable parallel programming with CUDA[J].Queue,2008,6(2):40-53.

[10] Kepler:The World’s Fastest,Most Efficient Hpc Architecture[EB/OL].[2012-03-16].http://www.nvidia.com/object/nvidia-kepler.html.

[11] 郭威.CT不完全投影數據重建算法研究[D].吉林大學,2011.

[12] Bell N, Garland M. Efficient sparse matrix-vector multiplication on CUDA[R].NVIDIA Technical Report NVR-2008-004,NVIDIA Corporation, 2008.

[13] Cook S. CUDA Programming: A Developer’s Guide to Parallel Computing with GPUs[M].Newnes,2012.

GPU ACCELERATED AND MATRIX OPTIMISED MEDICAL IMAGING RECONSTRUCTION

Zhou Wei1Dong Peng2Sun Yequan2Li Yaowu2

1(EquipmentDivision,WeifangMedicalUniversity,Weifang261053,Shandong,China)2(DepartmentofMedicalImaging,WeifangMedicalUniversity,Weifang261053,Shandong,China)

AbstractTo reconstruct the original image from huge projection data of CT (computed tomography) and other medical imaging equipment is a computationally intensive and time-consuming process. In order to accelerate CT image reconstruction, on the mainstream personal computer platform, we applied the CUDA (compute unified device architecture) technology, which is based on GPU (graphic processing unit) parallel computing, to an SIRT (simultaneous iterative reconstruction technique) algorithm ‘Cimmino’. Moreover, according to the characteristics of the GPU parallel architecture, we used CSR-ELL joint storage format to optimise data storage to further improve parallel efficiency. Experimental results showed that the GPU parallel Cimmino reconstruction method proposed in the paper could reduce 72% of the image reconstruction time, and improved the reconstruction speed under the premise of quality assurance.

KeywordsImage reconstructionCUDAComputed tomographySparse matrix

中圖分類號TP311.52

文獻標識碼A

DOI:10.3969/j.issn.1000-386x.2016.01.053

收稿日期：2014-05-25。山東省高等學校教學改革研究項目(2009074)；山東省研究生教育創新計劃項目(SDYC11085，SDYC13057)；中華醫學會醫學教育分會醫學教育研究項目(2012-SY-29)。周偉，助理實驗師，主研領域：并行計算。董鵬，教授。孫業全，教授。李耀武，實驗師。