一種利用Spark-GPU加速CT圖像重建的設(shè)計

熊威， 曾有靈， 李喆

(南方醫(yī)科大學(xué) 生物醫(yī)學(xué)工程學(xué)院， 廣東 廣州 510515)

CT(computed tomography, CT)圖像重建速度是衡量CT系統(tǒng)的重要指標(biāo).使用高性能計算硬件是加速圖像重建的重要手段，目前主要有使用CPU集群系統(tǒng)重建法、使用FPGA重建技術(shù)、使用GPU加速重建法以及使用細胞寬帶引擎(cell broadband engine,CBE)重建法[1].GPU擁有數(shù)以千計的計算核心，能同時創(chuàng)建成千上萬的線程，非常適合處理圖像重建中巨大的計算量.許多學(xué)者對GPU加速CT圖像重建進行了研究，Yan等[2]使用GPU加速FDK算法，并通過循環(huán)渲染到紋理、結(jié)合Z軸對稱和多個渲染目標(biāo)技術(shù)降低了切片幾何映射和投影的計算成本，提升了重建速度.Lu等[3]使用GPU加速同步代數(shù)迭代算法，采用射線驅(qū)動投影和體素驅(qū)動反投影技術(shù)使重建速度和常用的濾波反投影算法的速度相當(dāng).

統(tǒng)一設(shè)備計算框架(compute unified device architecture,CUDA)是英偉達(NVIDIA)公司推出的支持NVIDIA GPU的通用并行計算架構(gòu)，隨著GPU性能的提高和CUDA的成熟，GPU集群并行計算成為研究熱點，F(xiàn)an等[4]率先開發(fā)出可擴展的GPU集群，用于高性能科學(xué)計算和大規(guī)模仿真.消息傳遞接口MPI(message passing interface,MPI)是目前廣泛使用的并行計算開發(fā)環(huán)境.基于MPI-GPU[5-6]方案已經(jīng)成為設(shè)計GPU高性能計算集群的主流選擇.

大數(shù)據(jù)時代，傳統(tǒng)的并行計算框架已經(jīng)無法滿足快速、高效率的要求，Spark[7]正是在這種背景下誕生的產(chǎn)物.與CPU相比，GPU有眾多線程，在進行數(shù)據(jù)密集型計算時具有速度優(yōu)勢.Spark-GPU能夠同時使用CPU和GPU兩種計算資源，其性能高于以CPU為計算核心的Spark，因此利用GPU加速Spark應(yīng)用具有良好的應(yīng)用前景.

1 研究現(xiàn)狀

當(dāng)前兩大通用GPU計算框架分別是OpenCL(open computing language)和CUDA，實現(xiàn)Spark-GPU有兩種方案：Spark+OpenCL和Spark+CUDA.第一種方案具有代表性的是Segal等[8]提出的SparkCL框架，該框架包括Aparapi和Spark兩部分，其中由Aparapi執(zhí)行Java編寫的核函數(shù)并與OpenCL進行關(guān)聯(lián).與OpenCL相比，CUDA抽象程度高，易于使用，因此Spark-GPU多采用Spark+CUDA的方案.

在JVM(Java virtual machine)環(huán)境下，Spark通過JNI(Java native interface)與其他語言編寫的GPU核函數(shù)進行交互從而實現(xiàn)了對GPU的調(diào)用，如文獻[9-10].PyCUDA和Numba是CUDA支持Python所依賴的兩個主要的工具.有較多的研究者使用PyCUDA作為在Spark與GPU結(jié)合的主要工具，如文獻[11-12].

基于GPU加速的Spark框架能適合不同的計算任務(wù)，Li等[9]提出了針對機器學(xué)習(xí)的CPU-GPU異構(gòu)Spark計算平臺—HeteroSpark，Yuan等[10]利用GPU加速機器學(xué)習(xí)和SQL查詢，其中機器學(xué)習(xí)的速度提升了16.13倍，SQL查詢的速度提升了4.83倍.周情濤等[13]探討了Spark-GPU的實現(xiàn)，并闡述了利用PyCUDA實現(xiàn)K-means算法的流程.Ohno等[14]利用GPU 加速Spark中部分行動(action)和轉(zhuǎn)換(transform)操作.Serrano[11]在Spark框架中利用GPU加速重建高分辨率三維圖像(像素分辨率為2 048×2 048×2 048)，大幅縮短了重建時間.

利用GPU加速CT圖像重建是提高速度的重要方法，在高分辨率圖像重建[15]和批量圖像重建上都有良好的應(yīng)用.本文在此研究基礎(chǔ)上將Spark-GPU用于批量CT圖像重建，通過加速濾波反投影(fltered backprojection,FBP)算法和同時代數(shù)迭代重建技術(shù)(simultaneous algebraic reconstruction technique,SART)算法，使圖像的重建速度得到顯著提升.

2 Spark-GPU并行計算框架

2.1 Spark框架

Spark是加州大學(xué)伯克利分校的AMP實驗室(UC Berkeley AMP lab)所開源的通用并行計算框架，在大數(shù)據(jù)計算領(lǐng)域有著廣泛的應(yīng)用，如網(wǎng)絡(luò)流量分析，用戶行為預(yù)測，廣告推薦等業(yè)務(wù)上.Spark的通用性和易用性使其在生物醫(yī)學(xué)領(lǐng)域也有著良好的應(yīng)用前景，如人類基因組測序云平臺[16]、基因組學(xué)大數(shù)據(jù)分析[17]、醫(yī)學(xué)圖像并行處理[18]、批量醫(yī)學(xué)圖像重建[19]等.

Spark運行架構(gòu)包括運行在主節(jié)點(master)上的集群資源管理器(cluster manager)、運行作業(yè)任務(wù)的工作節(jié)點或從節(jié)點(worker node、slave node)、每個應(yīng)用的任務(wù)控制節(jié)點(driver)和每個工作節(jié)點上負責(zé)具體任務(wù)的執(zhí)行進程(executor).其中，集群資源管理器可以是Spark自帶的資源管理器，也可以是YARN(yet another resource negotiator)或Mesos等資源管理框架.Spark定義了兩大類方法：轉(zhuǎn)換(transform)和行動(action).轉(zhuǎn)換操作是惰性的，即通過構(gòu)造有向無環(huán)圖記錄彈性分布式數(shù)據(jù)集(resilient distributed data set,RDD)間的依賴關(guān)系，在行動操作時才會進行RDD變換.首先由SparkContext構(gòu)建有向無環(huán)圖(directed acyclic graph,DAG)，作業(yè)調(diào)度模塊(DAG scheduler)將DAG圖分解成階段(stage)，任務(wù)調(diào)度模塊(TaskScheduler)則將stage分解為task.任務(wù)控制節(jié)點向集群管理器申請資源，啟動Executor并向其發(fā)送應(yīng)用程序代碼和文件，由Executor執(zhí)行任務(wù)并將最終結(jié)果返回給任務(wù)控制節(jié)點.

2.2 Spark-GPU架構(gòu)

Spark-GPU總體架構(gòu)如圖 1所示，在Spark集群的基礎(chǔ)上，每個節(jié)點均結(jié)合了CPU和GPU兩種計算單元.當(dāng)一個節(jié)點內(nèi)有多個GPU時，為了充分發(fā)揮多GPU的性能，需要對GPU進行調(diào)度以選擇最合適的GPU進行計算.但在本架構(gòu)中，一個節(jié)點只配置一個GPU故無須調(diào)度.主節(jié)點將任務(wù)分發(fā)給各工作節(jié)點，各節(jié)點從本機硬盤讀取數(shù)據(jù)，在節(jié)點內(nèi)部由CPU負責(zé)復(fù)雜的串行邏輯計算，GPU負責(zé)數(shù)據(jù)密集的并行計算，通過兩者結(jié)合能充分提升集群性能.

圖1 Spark-GPU并行架構(gòu)

Spark-GPU程序可看作單機程序和Spark程序的結(jié)合，其中Spark程序負責(zé)讀取數(shù)據(jù)和保存數(shù)據(jù)，而單機程序是進行計算的部分.本文設(shè)計的程序分為3個部分：①主程序使用thunder讀取數(shù)據(jù)并創(chuàng)建RDD， RDD中的每個元素為一個完整的投影數(shù)據(jù).在創(chuàng)建RDD時可以設(shè)置n個元素為1個task，默認task數(shù)等于元素數(shù).當(dāng)n較大時，執(zhí)行一個task會消耗更多內(nèi)存.主程序還負責(zé)將各節(jié)點返回的結(jié)果保存到本地. ②圖像重建程序負責(zé)接收RDD中的每個元素并進行重建，通過網(wǎng)絡(luò)將重建結(jié)果返回給Master. ③圖像重建程序?qū)⒉糠謹(jǐn)?shù)據(jù)傳到GPU顯存，調(diào)用GPU加速接口執(zhí)行計算任務(wù).主節(jié)點負責(zé)第一部分并不參與計算，但仍需要配置顯卡，從節(jié)點負責(zé)執(zhí)行①和②.

2.3 開源工具Numba

早期CUDA支持C、C++和Fortran語言，隨著Python在科學(xué)計算領(lǐng)域的突出表現(xiàn)，NVIDIA公司先后推出PyCUDA[20]和Numba[21]兩款工具，使Python成為CUDA支持的第4種語言.PyCUDA是支持RTCC(GPU run-time code generation)技術(shù)的Python開源工具.其原理是將C語言編寫的核函數(shù)源碼當(dāng)作字符串傳給Source Module的構(gòu)造函數(shù)，通過動態(tài)編譯生成機器代碼.PyCUDA采用了編譯器緩存機制，在重復(fù)調(diào)用核函數(shù)時能減少編譯次數(shù)，節(jié)省編譯時間.Numba是基于底層虛擬機(low level virtual machine，LLVM)的即時編譯工具.Numba提供和CUDA C類似的編程模型，稱為Numba CUDA模型，該模型包含了CUDA C的大部分功能，并且完全使用Python語言，支持在核函數(shù)內(nèi)部使用Cmath、 Math函數(shù)庫，因此Numba能充分發(fā)揮Python的靈活性和擴展性.Spark并不支持直接調(diào)用GPU，通過導(dǎo)入Numba工具庫可以直接用Python編寫核函數(shù)，從而在Spark中使用GPU進行計算.

3 利用GPU加速CT圖像重建分析

3.1 FBP 算法及并行化

CT重建算法按照其重建方式可分為兩類：解析重建算法和迭代重建算法.FBP是經(jīng)典的解析重建算法.根據(jù)中心切片定理，二維圖像可表達為：

t=xcosθ+ysinθ

(1)

f(x,y)為x-y平面上定義的密度函數(shù)，P(ω,θ)是密度函數(shù)的二維傅里葉變換與ωy軸夾角為θ的切片,|ω|是斜坡濾波器.

FBP算法的步驟可分為濾波和反投影兩步.濾波的時間復(fù)雜度為O(N2log2N),反投影的時間復(fù)雜度為O(N3)，可見反投影計算量大、耗時長.GPU濾波一般采用cuFFT庫，時間復(fù)雜度為O(Nlog2N).對短序列進行快速傅里葉變換(fast Fourier transform,FFT)FFT變換，GPU的速度要慢于CPU的速度.

GPU可將所有角度的反投影一步完成，時間復(fù)雜度減少為O(1)，可見將反投影并行化是加速FBP重建的關(guān)鍵.主機將每個角度濾波后的投影數(shù)據(jù)傳輸?shù)紾PU全局內(nèi)存，由GPU進行反投影并將各個角度的投影數(shù)據(jù)疊加，再將最終數(shù)據(jù)傳輸?shù)街鳈C內(nèi)存.例如180個角度的投影數(shù)據(jù)，每個角度的數(shù)據(jù)長度為725，將投影數(shù)據(jù)進行FFT變換和反投影的時間對比如圖 2.

圖2 FFT和反投影時間對比

由圖2可以看出在CPU端進行FFT的時間遠小于在GPU端進行FFT的時間，由于GPU內(nèi)核在首次啟動時耗時更長，因此在CPU端進行濾波是最優(yōu)的選擇.反投影計算量大，因此GPU多線程更有優(yōu)勢，利用GPU進行反投影的時間約為CPU的1/3.通過分析可知在主機端進行濾波，在GPU端進行反投影是最佳的組合.

3.2 SART算法及并行化

SART算法是將原圖像看作待求解的矩陣，并通過方程來求解.設(shè)原圖像大小是像素，如式(2)所示：

WX=P

(2)

X即為待求的N個像素值(N=n2)，是一維列向量.P是所有射線投影數(shù)據(jù)向量，P=(p1,p2,…,pM),M是投影射線數(shù)量.W是M×N維矩陣，其元素wij表示第j個體素對第i個投影數(shù)據(jù)pi的貢獻值.

同時代數(shù)迭代算法(SART)具體步驟如下：

1)對未知圖像進行賦初值：

(3)

2)計算圖像的投影值：

(4)

3)計算誤差：

(5)

4)計算修正值：

(6)

5)對第j個像素值進行修正：

(7)

λ為松弛因子.

6)以上迭代結(jié)果作初始值，重復(fù)2)-5)步驟.

SART可以簡化為4個步驟，即正投影、修正、反投影、更新，基于GPU加速的SART算法多采用體素驅(qū)動正投影，并行遍歷以及射線驅(qū)動反投影[22].正投影和反投影總時間復(fù)雜度均為O(N3)，并行化后時間總的復(fù)雜度變?yōu)?O(N)，復(fù)雜度降低，速度提升大.FBP中所有角度的反投影是同時進行的，而SART的反投影是逐次進行，一次子迭代過程包含4個步驟，只對一個角度的數(shù)據(jù)進行反投影，算法的差異導(dǎo)致SART的速度遠慢于FBP的速度.

根據(jù)兩種圖像重建算法共同的特點，本文設(shè)計了兩種利用GPU加速的函數(shù)：正投影函數(shù)和反投影函數(shù).將數(shù)據(jù)密集型計算任務(wù)轉(zhuǎn)移到GPU，能夠充分利用了GPU多線程，大幅提升程序運行速度.針對SART迭代計算的特點，使用GPU進行修正和更新，通過將部分?jǐn)?shù)據(jù)常駐GPU內(nèi)存，在調(diào)用核函數(shù)時能減少主機和GPU的數(shù)據(jù)傳輸次數(shù)，加快程序運行速度.實驗表明使用GPU緩存能夠使SART程序運行時間減少32%.

4 實驗

4.1 實驗環(huán)境配置

集群硬件配置如表 1所示.集群由一臺PC作為Master節(jié)點，3臺PC作為Slave節(jié)點，各節(jié)點通過百兆交換機連接.每個節(jié)點均運行Centos 7.0系統(tǒng)，系統(tǒng)所需環(huán)境為Spark-2.3，Numba-36.0，thunder-1.4.2，CUDA-9.1.

表1 集群硬件配置Table 1 Cluster hardware configuration

本文所用的CT圖像數(shù)據(jù)為Shepp-Logan模型，圖像像素分辨率分別為 256×256、512×512.投影角度0°～180°，共180個角度，每個角度投影下的射線數(shù)分別為363和725.SART算法松弛因子為0.02，迭代次數(shù)為2×180次.重建圖像效果如圖3.

圖3 原圖像、FBP重建圖像和SART重建圖像

4.2 實驗結(jié)果

本實驗一次性重建圖像數(shù)量分別為30，60，90，300，600，900，1 200張，圖像數(shù)量為3的整數(shù)倍，便于主節(jié)點分配.記錄3個節(jié)點下重建上述數(shù)量圖像的時間，詳細數(shù)據(jù)如表 2所示，重建時間以秒為單位，數(shù)據(jù)保留小數(shù)點后兩位.

表2 重建圖像所需的時間Table 2 Time required to reconstruct the image

使用SART和FBP算法重建一張圖像的時間隨數(shù)量的變化如圖4.可以看出，隨著圖像數(shù)量的增加，平均重建一張圖像所需的時間逐漸減少，速度越來越快.這是因為從提交Spark應(yīng)用開始，任務(wù)控制節(jié)點向集群管理器申請資源，任務(wù)調(diào)度，啟動executor以及集群間的通信等都需要一定的時間，這些時間都是系統(tǒng)開銷.在圖像數(shù)量較少時，系統(tǒng)開銷在總時間中占比大，隨著任務(wù)執(zhí)行時間變長，系統(tǒng)開銷相對減少，重建速度逐漸提升并趨于穩(wěn)定.FBP算法的時間最快為0.03 s，SART算法的時間最快為0.85 s.同樣參數(shù)下，在單機上用CPU運行FBP算法的時間為1.31秒，運行SART算法的時間為8.34 s，經(jīng)過Spark-GPU加速后，F(xiàn)BP算法的速度有近40倍的提升，SART算法的速度有近10倍的提升.由以上分析可知，隨著數(shù)據(jù)增多，Spark-GPU加速效果明顯，速度逐漸達到理想狀態(tài)，適合大規(guī)模計算任務(wù).

4.3 集群性能分析

通過計算相對加速比[23]可以判斷集群的性能和效率，其定義式為：

相對加速比反映了當(dāng)集群節(jié)點數(shù)目增加時集群性能提升的大小，在理想情況下，加速比與節(jié)點數(shù)量成正比.在重建900張分辨率為512×512的圖像時，集群的加速比如圖5所示，可以看出加速比與節(jié)點數(shù)量并不嚴(yán)格成正比，這是由于各節(jié)點的CPU性能不完全相同導(dǎo)致加速比降低.SART算法的運行時間遠高于FBP算法，系統(tǒng)開銷相對減少，所以SART算法的加速比大于FBP算法的加速比.使用性能相同的處理器、保證硬件配置一致，防止集群出現(xiàn)性能短板，合理地分配資源和設(shè)置任務(wù)大小等方法都可以提升加速比.

圖5 集群加速比

以Slave03節(jié)點為例，圖 6記錄該節(jié)點運行的前200 s內(nèi)CPU和內(nèi)存利用率，可以看出在程序開始后CPU的使用率接近100%，內(nèi)存使用率在30%上下浮動.在默認情況下，Spark給每個節(jié)點分配的任務(wù)數(shù)與該節(jié)點CPU的核數(shù)相同，這使得CPU性能被充分利用.例如某節(jié)點的CPU為四核，Spark可以發(fā)起4個線程同時執(zhí)行4個Task.Spark默認的任務(wù)調(diào)度策略為FIFO(first in first out)，一個task所需的內(nèi)存大小與該task的分區(qū)大小有關(guān)，當(dāng)一個task完成后，結(jié)果會立即返回給主節(jié)點并釋放內(nèi)存，從而避免內(nèi)存占用過高.

與從節(jié)點相比，由于主節(jié)點并不參與計算，而是接收各個節(jié)點傳回的數(shù)據(jù)，因此對網(wǎng)絡(luò)帶寬的要求高，圖 7記錄了FBP和SART程序開始后的45 s內(nèi)主節(jié)點網(wǎng)絡(luò)下行速度，可以看出FBP程序基本將網(wǎng)絡(luò)帶寬全部占用，這是因為FBP算法的速度更快，短時間內(nèi)大量的數(shù)據(jù)要傳回主節(jié)點.SART的速度較慢，從節(jié)點有充分的時間傳輸數(shù)據(jù)，因此主節(jié)點的下載速度較低，只出現(xiàn)有規(guī)律的波動.網(wǎng)絡(luò)帶寬對Spark任務(wù)調(diào)度亦有較大影響，由于FIFO策略的特點，后續(xù)任務(wù)等待被執(zhí)行直至從節(jié)點將結(jié)果傳給主節(jié)點，當(dāng)網(wǎng)絡(luò)帶寬較小時，任務(wù)調(diào)度時間會明顯變長.在FBP程序中，任務(wù)調(diào)度時間和執(zhí)行時間相當(dāng)導(dǎo)致計算效率嚴(yán)重降低.通過提升帶寬能縮短傳輸時間，或者將數(shù)據(jù)直接保存在從節(jié)點的磁盤中，后者能將總時間縮短50%.

主節(jié)點接收傳回的數(shù)據(jù)，因此內(nèi)存占用大于從節(jié)點.Spark對內(nèi)存的依賴使得內(nèi)存大小直接影響了能夠處理數(shù)據(jù)的多少，為了能處理更多的數(shù)據(jù)，可以將數(shù)據(jù)集拆分，分批創(chuàng)建RDD，待前一個RDD任務(wù)完成后再執(zhí)行下一個RDD任務(wù).除了集群硬件性能外，Spark系統(tǒng)配置參數(shù)的設(shè)置也會對任務(wù)的運行性能有很大的影響[24]，如maxResultSize太小會導(dǎo)致driver沒有足夠空間容納序列化結(jié)果，太大會導(dǎo)致driver內(nèi)存溢出.

圖6 slave03系統(tǒng)資源使用情況

圖7 Master節(jié)點下行速度

5 討論

大規(guī)模圖像重建對計算機性能有了更高的要求，本文構(gòu)建由4節(jié)點組成的Spark-GPU集群，實現(xiàn)CPU和GPU結(jié)合的并行計算，大幅提升了圖像重建的速度，并可通過增加節(jié)點的方式獲得更大的加速比以滿足更高的任務(wù)要求.批量圖像重建是將數(shù)據(jù)由大化小，分批重建，其核心思想是MapReduce，在超高分辨率圖像和三維圖像的加速重建中也有應(yīng)用價值.本文使用Spark-CUDA-Numba-Thunder的結(jié)合模式是為整合Spark與GPU提供了新方向.

由于各節(jié)點通過百兆交換機相連，不能滿足大數(shù)據(jù)量的傳輸要求，網(wǎng)絡(luò)配置還需進一步升級，同時設(shè)計新的調(diào)度策略提升數(shù)據(jù)的傳輸效率，降低調(diào)度時間.GPU內(nèi)存一般不及主機內(nèi)存大且不可擴展，如何優(yōu)化算法減少GPU內(nèi)存使用，對充分發(fā)揮GPU性能具有重要意義.Spark還無法在應(yīng)用間管理和調(diào)度GPU資源，在應(yīng)用內(nèi)對GPU的使用和單機程序并無區(qū)別，當(dāng)一個節(jié)點內(nèi)有多個GPU時，對GPU的管理調(diào)度十分重要.有文獻提出將GPU作為單獨的計算資源進行調(diào)度，無論在應(yīng)用間和應(yīng)用內(nèi)都能自動選擇最合適的GPU[10-11]，但需要額外設(shè)計GPU的調(diào)度程序.如何將GPU融合到Spark框架中，不需要外部工具包并在多語言平臺實現(xiàn)統(tǒng)一是重要的研究方向.