基于GPU加速和矩陣優(yōu)化的醫(yī)學圖像重建

周　偉　董　鵬　孫業(yè)全　李耀武

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

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

?

基于GPU加速和矩陣優(yōu)化的醫(yī)學圖像重建

周偉1董鵬2孫業(yè)全2李耀武2

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

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

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

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

0引言

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

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

1圖像重建原理

1.1圖像重建數(shù)據(jù)采集

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

Ax=b

(1)

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

圖1　系統(tǒng)矩陣示意圖

1.2Cimmino重建算法

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

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

(2)

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

(3)

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

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

2.1串行Cimmino算法

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

算法1串行Cimmino算法

(1) for k ← 0:K

//K為最大迭代次數(shù)

(2) M_dp_rxk ← M(b - Axk)

(3) xk+1← xk+ λATM_dp_rxk

(4) end

2.2并行的Cimmino算法

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

算法2并行Cimmino算法

(1) for k ← 0:K

//K為最大迭代次數(shù)

(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存儲結構優(yōu)化

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

圖2　系統(tǒng)矩陣CSR-ELL存儲優(yōu)化

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

3實驗方法和實驗結果

3.1實驗方法

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

3.2實驗結果

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

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

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

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

4結語

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

參考文獻

[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不完全投影數(shù)據(jù)重建算法研究[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)；山東省研究生教育創(chuàng)新計劃項目(SDYC11085，SDYC13057)；中華醫(yī)學會醫(yī)學教育分會醫(yī)學教育研究項目(2012-SY-29)。周偉，助理實驗師，主研領域：并行計算。董鵬，教授。孫業(yè)全，教授。李耀武，實驗師。